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ABSTRACT 

This  report  details  the  areas  in  which  signal  processing  techniques  have  enabled  DSTO 
Sydney  to  enhance  its  existing  magnetic  testing  capability,  and  present  viable  data  processing 
options  for  their  implementation.  This  entails  not  only  accurate  sensor  calibration,  but  also 
calibration  of  the  excitation  coils  used  within  the  magnetic  test  system.  Reduction  of  external 
noise  influences  through  the  use  of  high  permeability  shielding  and  filtering  are  detailed, 
presenting  modelling  of  coil-shield  interactions,  and  experimental  tests  for  various  filter 
characteristics.  Various  magnetic  coil  designs,  including  DSTO's  previously  undocumented 
"Ternan"  coil,  are  also  investigated,  leading  to  the  design  of  a  new  magnetic  coil,  the 
"ELFcage"  coil,  being  developed.  Field  uniformity  measurements  for  the  "Ternan"  magnetic 
coil  are  also  presented. 
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Magnetic  Test  Facility  -  Sensor  and  Coil  Calibrations 


Executive  Summary 


This  report  contains  the  author's  Master  of  Sciences  (Defence  Signals  and  Information 
Processing)  thesis,  completed  as  part  of  the  DSTO  Continuing  Education  Initiative. 

This  thesis  investigates  the  coil  geometries  of  the  Australian  Defence  Science  and 
Technology  Organisation's  previously  undocumented  "Ternan"  magnetic  coil  facility. 
Computational  modelling  has  been  conducted  to  investigate  the  extent  of  the  usable 
volume  of  uniform  magnetic  field  generated  by  the  magnetic  coils,  as  well  as 
modelling  of  an  enhanced  coil  design  named  "ELFcage"'.  Experimental  measurements 
have  found  the  Ternan  coil  to  have  a  very  large  area  of  uniformity,  permitting  the 
testing  of  objects  previously  deemed  too  large  to  be  easily  investigated. 

Also,  this  paper  investigates  the  impact  high  permeability  magnetic  shielding  material 
has  on  an  enclosed  magnetic  coil.  Computational  modelling  has  shown  the  shielding 
material  enhanced  the  magnitude  and  the  uniformity  of  the  magnetic  flux  density 
within  the  coil. 

Finally  this  paper  investigates  methods  used  for  calibrating  magnetic  sensors  and 
triaxial  magnetic  coil  systems,  to  enable  the  development  of  an  automated  calibration 
procedure  for  triaxial  magnetic  coils.  Since  existing  magnetic  coil  calibration  techniques 
found  in  literature  do  not  fully  account  for  secondary  errors  such  as  orthogonality,  the 
TWOSTEP  magnetometer  calibration  algorithm  is  recast  to  enable  its  use  to  calibrate 
triaxial  magnetic  coils.  Experimental  measurements  have  found  the  procedure  enables 
an  accurate  coil  calibration  to  be  conducted  with  90  s  of  measurement  data, 
significantly  reducing  the  error  in  the  generated  magnetic  field  when  using 
theoretically  derived  calibration  parameters 

Given  the  DSTO  magnetic  test  facility  is  already  built  and  available  for  use  it  is 
recommended  that  the  existing  Ternan  coil  configuration  be  kept,  and  be  used  for 
characterisation  of  ranging  system  sensors,  and  measurement  of  magnetic  signatures  of 
bulky  equipment.  It  is  also  recommended  a  small  2-axis  ELFcage  coil  be  built  for  y-axis 
and  z-axis  field  stimulation  of  small  sensors  within  a  shielded  laboratory  environment, 
for  use  in  the  design  and  development  of  future  ranging  systems,  and  that  an 
additional  removable  solenoid  be  built  to  stimulate  the  x-axis  field. 
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Abstract 


This  thesis  investigated  the  methods  used  for  calibrating  magnetic  vector  sensors  and 
triaxial  magnetic  coils,  to  enable  the  development  of  an  automated  calibration  pro¬ 
cedure,  as  well  as  investigating  the  coil  geometries  of  the  previously  undocumented 
"Ternan"  coil,  as  part  of  its  ongoing  refurbishment. 

Since  existing  magnetic  coil  calibration  techniques  found  in  literature  did  not  fully 
account  for  secondary  errors  such  as  orthogonality,  a  magnetometer  calibration  algo¬ 
rithm  was  extended  to  enable  its  use  to  calibrate  triaxial  magnetic  coils. 

Computational  modelling  was  also  conducted  to  investigate  the  extent  of  the  usable 
volume  of  uniform  magnetic  field  generated  by  the  magnetic  coils,  and  the  impact  high 
permeability  materials  would  have  on  it,  if  used  to  shield  the  system  from  external 
magnetic  influences. 

The  extended  calibration  algorithm  was  found  to  be  effective  at  determining  the  mag¬ 
netic  coil  fabrication  errors,  allowing  them  to  be  counteracted  by  the  use  of  a  calibration 
matrix  within  the  control  software  of  the  system. 

The  Ternan  coil  was  found  to  have  a  very  large  area  of  uniformity,  permitting  the  test¬ 
ing  of  objects  previously  deemed  too  large  to  be  easily  investigated;  and  the  volume  of 
magnetic  uniformity  was  found  to  be  extendable  if  appropriate  shielding  was  placed 
around  the  system. 
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Introduction 


Mines /Torpedoes1  were  first  used  against  ships  in  late  Dec  1777  -  early  Jan  1778,  in 
the  Delaware  River  at  Philadelphia;  kegs  filled  with  gunpowder  which  were  set  to 
explode  when  a  flintlock  trigger  fired  on  contact  with  targets  that  drifted  into  them 
once  they  were  released  down  a  river  (Hartmann  1979).  Since  then  mines  progressively 
evolved  to  include  more  complicated  triggering  mechanisms,  providing  both  greater 
safety  during  deployment  and  better  target  discernment.  Contact  mines  such  as  these 
were  used  throughout  the  First  World  War,  and  more  recently  by  countries  such  as 
Iran. 

The  Second  World  War  saw  the  development  of  magnetic,  pressure,  and  acoustic  mines. 
Their  evolution  has  led  to  the  development  of  mines  capable  of  detecting  and  classify¬ 
ing  vessels  based  on  increasingly  more  detailed  characteristics  of  their  signatures,  i.e., 
the  physical  characteristics  of  the  platform  for  the  purpose  of  detection  and  classifica¬ 
tion.  Consequently  military  organisations  and  scientists  globally  were  and  are  endeav¬ 
oring  to  minimise  signatures  as  far  as  possible,  whilst  concurrently  researching  other 
methods  to  counter  the  continually  changing  mine  threat.  The  determination  of  these 
signatures  is  called  "ranging";  and  this  can  be  performed  on  a  variety  of  platforms, 
varying  from  large  vessels  down  to  smaller  platforms  such  as  autonomous  underwa¬ 
ter  vehicles  (AUVs)  or  diving  gear. 

The  drive  to  accurately  measure  vessel  signatures  has  led  to  the  use  of  more  accurate 
sensors  able  to  detect  absolute  pressure  signals  to  part-per-billion  resolutions  (Schaad 
2009).  Acoustic  sensors  have  progressed  to  the  point  where  they  are  limited  by  envi¬ 
ronmental  noise  (Stansfield  1991)  and  magnetic  sensors  are  capable  of  detecting  changes 

1  Mines  were  initially  called  Torpedoes 
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in  the  sub  nano-Tesla  (nT)  ranges  (Lenz  and  Edelstein  2006).  Because  these  resolutions 
are  so  small,  even  slight  errors  in  magnetometer  manufacture,  such  as  a  misalignment 
of  0.1  degrees,  could  lead  to  several  tens  of  nT  errors  in  the  readings,  which  can  easily 
swamp  the  signals  being  measured. 

The  three  main  sensors  used  in  ranging  are  acoustic,  pressure  and  magnetic.  DSTO 
MOD  currently  has  the  ability  to  test  acoustic  sensors  at  Woronora  Dam  (Sydney)  and 
pressure  sensors  in-house.  However  magnetic  sensors  have  not  been  able  to  be  prop¬ 
erly  characterised  since  the  shutdown  of  the  Maritime  Operations  Division  of  DSTO 
Maribyrnong  (Melbourne)  began  in  2004,  during  which  time  the  magnetic  testing  and 
calibration  equipment  was  relocated  into  storage  in  Sydney.  In  this  project  the  pre¬ 
existing  magnetic  coil  system  from  Maribyrnong  (referred  to  as  the  Ternan  coil  mag¬ 
netic  volume,  shown  in  Figure  1.1)  has  since  been  thoroughly  modelled  to  gain  a  com¬ 
prehensive  understanding  of  the  coil  and  magnetic-field  topology,  prior  to  being  rein¬ 
stated  as  a  fully-operational  magnetic-field  calibration  facility.  Furthermore  an  alterna¬ 
tive  magnetic  coil  design  (referred  to  as  the  EFFcage)  was  also  thoroughly  investigated 
and  was  shown  to  exhibit  a  higher  degree  of  magnetic  field  uniformity  compared  to 
the  Ternan  magnetic  volume.  This  was  a  significant  finding. 


Figure  1.1:  Ternan  coil  magnetic  volume  as  delivered  to  Sydney  in  2004  (Coil  length 
=  8  m,  coil  diameter  =  2  m). 
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This  thesis  project  will  present  the  signal  processing  techniques  which  are  being  used 
in  the  re-establishment  of  a  large-volume  magnetic  testing  and  calibration  facility  using 
the  Ternan  coil,  and  their  implementation  for  use  in  a  small  scale  coil  located  within  a 
magnetically  noisy  laboratory 

The  thesis  will  commence  with  an  overview  of  the  types  of  sensors  available,  their 
sources  of  error,  and  calibration  techniques  used  to  calibrate  them.  An  overview  of 
magnetic  coils  will  then  introduce  the  advantages  of  using  the  non-standard  coil  de¬ 
sign  utilised  in  the  Ternan  coil  system,  as  well  as  the  expected  sources  of  errors  and 
limitations  of  calibration  techniques  available  within  the  literature.  Noise  sources  and 
the  effectiveness  of  magnetic  shielding  will  then  be  addressed,  before  assessing  suit¬ 
ability  of  two  estimators  (the  "TWOSTEP"  and  "Geometric  Approach"  methods)  and 
their  performance.  In  this  context  the  TWOSTEP  formulation  was  re-derived  to  ver¬ 
ify  the  equations  and  was  coded  in  MATLAB,  enabling  processing  performance  to  be 
tested  under  various  noise  levels  and  sample-set  sizes.  Finally  a  series  of  field  mea¬ 
surements  will  be  presented  for  both  sensors  and  coils. 
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2.1  Sensors 


Magnetic  sensors  can  be  divided  into  two  classes,  scalar  magnetometers  and  vector 
magnetometers. 

As  the  name  indicates,  scalar  magnetometers  measure  the  total  field  strength  at  a  point, 
whilst  vector  magnetometers  provide  the  3D  components  making  up  the  field.  Vector 
magnetometers  generally  have  higher  sensitivity  and  enable  higher  sampling  frequen¬ 
cies  than  their  scalar  counterparts,  with  the  exception  of  SQUIDs  (Superconducting 
Quantum  Interference  Devices),  which  are  bound  to  low  sampling  frequencies.  Vec¬ 
tor  magnetometers  are  also  smaller  than  the  scalar  sensor,  since  they  do  not  require 
volumes  of  liquids  or  gases  to  excite. 

Scalar  magnetometers  are  less  sensitive,  bulkier  and  provide  lower  sampling  frequen¬ 
cies  than  vector  magnetometers,  but  their  use  becomes  apparent  once  fabrication  er¬ 
rors  (e.g.  sensor  orthogonality  and  bias  voltage)  within  vector  sensors  are  taken  into 
consideration. 

Selection  of  the  magnetic  sensor  for  a  system  will  often  have  a  major  impact  on  the 
system's  final  errors  and  capabilities;  as  such  it  is  important  to  identify  the  most  ap¬ 
propriate  sensor  for  its  intended  purpose(s). 

Lenz  and  Edelstein  (2006)  and  Tumanski  (2007)  provide  good  overviews  of  existing 
and  emerging  vector  and  scalar  magnetic  sensors,  a  few  of  which  are  summarised  in 
Sections  2.1.1  and  2.1.2. 
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2.1.1  Vector  magnetometers 

2. 1.1.1  Search-coil  magnetometers 

Search-coils  consist  of  a  conducting  wire  wound  around  a  core  commonly  made  of 
high  permeability  materials.  The  response  is  governed  by  Faraday's  law  of  induction, 
where  a  voltage  is  generated  proportional  to  the  change  of  flux  through  a  coil.  Signals 
in  the  tens  of  FT  have  been  detected  with  these  sensors,  but  search-coils  can  also  be 
designed  to  detect  extremely  strong  signals,  which  are  orders  of  magnitude  stronger 
(i.e.,  >  1  mT).  The  typical  frequency  of  operation  is  1  Hz  to  1  MHz,  dependent  on  the 
inductance,  stray  capacitance  and  resistance  of  the  system.  These  sensors  are  unlikely 
to  be  useful  in  the  measurement  of  sub  1  Hz  signals  due  to  the  small  amplitude  of 
signals  induced  into  the  sensor. 

2. 1.1. 2  Fluxgate  magnetometers 

Fluxgate  sensors  consist  of  a  ferromagnetic  core,  around  which  drive  and  sense  coils 
are  wound.  They  use  magnetic  induction  coupled  with  the  change  of  permeability  of 
the  core  as  it  goes  into  saturation  at  high  fields  caused  by  the  applied  drive  signal,  to 
induce  even-numbered  harmonics  of  the  drive  signal  into  the  sense  coil.  The  ampli¬ 
tude  of  the  second  harmonic  is  proportional  to  the  ambient  field  strength  in  the  axial 
direction  of  the  coil.  Sensitivity  range  is  from  ~  10  pT  to  ~  10  mT.  The  sampling 
frequency  of  these  sensors  is  limited  by  the  driving  signal,  which  is  in  turn  limited  to 
~  10  kHz  by  the  response  time  of  the  core,  thereby  enabling  flux  variations  in  the  kHz 
range  to  be  measured.  Fluxgate  magnetometers  do  not  have  a  lower  limit  to  their  fre¬ 
quency  response,  and  hence  can  be  used  to  measure  frequencies  all  the  way  down  to 
DC  fields. 

2. 1.1. 3  SQUID  magnetometers 

The  SQUID  sensor  consists  of  a  small  superconducting  loop  in  which  a  quantized  cur¬ 
rent  passes,  dependent  on  the  flux  density  passing  through  the  loop.  Superconduc¬ 
tivity  is  achieved  at  low  temperatures,  requiring  liquid  Helium  —  referred  to  as  low 
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temperature  superconductivity  (<  4  K),  or  liquid  nitrogen  —  referred  to  as  high  tem¬ 
perature  superconductivity  (77  K),  dependent  on  the  superconductive  material  used. 
The  current  is  measured  by  observing  the  pattern  of  the  oscillating  current  passing 
through  a  "weak  link"2,  allowing  information  to  be  inferred,  due  to  its  similarity  to 
interference  patterns  in  light,  but  in  this  case  related  to  flux  quantization.  Sensitivity 
is  limited  by  magnetic  field  noise  to  ~  10  fT.  The  devices  can  easily  be  set  to  null  out 
static  fields,  making  them  well  suited  for  anomaly  detection.  Frequency  response  is 
from  DC  to  ~  1Hz.  SQUIDs  can  be  configured  to  act  as  gradiometers,  but  at  the  ex¬ 
pense  of  resolution.  The  requirement  of  a  supply  of  liquid  Nitrogen  or  Helium  is  a 
major  impediment  to  their  widespread  use  in  the  field. 


2. 1.1. 4  Hall  Effect  magnetometers 

These  sensors  work  on  the  Hall  effect,  a  property  of  conductors  and  semiconductors  by 
which  a  strong  magnetic  field  perpendicular  to  the  plane  of  a  rectangle  of  conductive 
material  with  current  passing  along  its  length  will  generate,  due  to  the  Lorenz  force,  a 
voltage  across  the  width  of  the  rectangle,  proportional  to  the  field  strength.  Sensitivity 
typically  ranges  from  ~  100  nT  to  ~  0.1  T,  and  can  sample  signals  up  to  1  MHz. 


2. 1.1. 5  Other  vector  magnetometers 

A  variety  of  other  vector  magnetometers  are  available,  including  MEMS  (Micro  Elec¬ 
tro  Mechanical  Systems)  devices,  magneto-resistive  devices  (which  exploit  the  change 
in  resistance  of  certain  materials  in  a  magnetic  field),  magneto  optical  sensors  (which 
rotate  the  plane  of  a  polarised  light  passing  through  a  field),  and  magneto-strictive  sen¬ 
sors  (which  use  the  physical  change  of  size  of  materials  as  the  applied  field  changes). 

Most  vector  magnetometers  exhibit 1  //  noise  and  are  susceptible  to  errors  introduced 
by  rotational  vibrations. 


2a  thin  layer  of  insulator,  or  a  thinning  of  the  superconducting  loop 
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2.1.2  Scalar  magnetometers 

Scalar  sensors  are  insensitive  to  rotational  vibrations,  and  will  reliably  measure  the  to¬ 
tal  field  strength  at  a  point  but,  due  to  the  physical  processes  involved  in  these  devices, 
they  are  only  usable  for  sampling  frequencies  below  ~  10Hz,  and  require  the  total  field 
strength  to  be  above  a  minimum  level. 

2. 1.2.1  Optically  pumped  magnetometers 

These  sensors  detect  the  splitting  of  spectral  lines  when  Cesium  and  Rubidium  atoms, 
for  example,  are  placed  in  a  magnetic  field.  Detection  is  achieved  by  tracking  the  fre¬ 
quency  of  a  variable  radio  frequency  (RF)  source,  which  is  used  to  release  excited  elec¬ 
trons  from  their  high  energy  state;  this  is  achieved  via  a  feedback  loop  involving  a 
photo-sensor  to  keep  the  RF  source  tuned  to  produce  maximum  light  intensity.  Sen¬ 
sitivity  typically  ranges  from  1  pT  to  100  y.T.  These  devices  exhibit  "dead  zones", 
requiring  multiple  sensors  to  obtain  a  full  omnidirectional  spherical  response. 

2. 1.2. 2  Nuclear  precession  magnetometers 

Also  known  as  proton  precession  magnetometers,  nuclear  precession  magnetometers 
temporarily  align  protons  in  a  hydrocarbon  fluid  to  a  magnetic  field  generated  by  a 
coil.  When  the  induced  field  is  removed,  the  protons  precess  about  the  ambient  field, 
creating  a  signal  in  the  coil,  the  frequency  of  which  is  dependent  on  the  ambient  field 
strength.  Sensitivity  typically  ranges  from  0.1  nT  to  100  }iT.  Nuclear  precession  mag¬ 
netometers  have  no  dead  zones. 

2. 1.2. 3  Overhauser  magnetometers 

These  sensors  are  an  evolution  of  the  nuclear  precession  magnetometers,  which  use  a 
coupling  between  the  electron  spins  and  the  protons,  to  get  a  factor  of  1000  increase  in 
the  nuclear  polarisation,  hence  providing  better  noise  levels,  as  low  as  15  PT/ vTtz  at  1 
Hz. 
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2.1.3  Sources  of  error 

Small  instrumentation  errors  are  generally  introduced  into  the  magneto-sensing  device 
by  the  fabrication  process  itself.  For  vector  magnetometers  these  errors  include  sensor 
bias  /offset,  sensor  drift,  sensor  orthogonality  and  coordinate  system  misalignment. 
Obviously  a  scalar  magnetometer  cannot  experience  orthogonality  or  coordinate  sys¬ 
tem  errors,  since  it  does  not  contain  coordinate  information,  but  as  mentioned  earlier 
it  can  suffer  from  other  errors,  including  the  "dead  zone",  due  to  which  the  sensitivity 
response  of  the  sensor  is  not  a  spherical  beam  pattern. 

Reduction  of  many  of  these  errors  can  be  achieved  through  sensor  calibration,  a  pro¬ 
cess  which  initially  involved  direct  measurement  of  offset  and  the  linear  approxima¬ 
tion  of  sensitivities,  but  which  has  matured  to  utilise  a  swathe  of  signal  processing 
techniques  including  Kalman  filters,  neural  networks,  maximum  likelihood  estimators 
and  least  mean  square  estimators,  as  well  as  image  processing  techniques  such  as  el¬ 
lipsoid  fitting.  What  follows  is  an  overview  of  techniques  used  in  recent  papers,  the 
most  relevant  of  which  will  be  further  investigated  for  implementation  in  the  Ternan  & 
ELFcage  coil  systems.  This  technique  can  be  extended  for  use  with  other  coil  designs. 


2.1.4  Sensor  calibration 

IEEE  standards  propose  several  ways  of  calibrating  a  sensor  (Frix  et  al.  1994,  Mirza- 
eva  et  al.  2012),  namely  using  a  single  square  coil  as  a  reference  source,  using  a  Helmholtz 
coil  as  a  reference  source,  and  calibrating  against  a  reference  sensor.  While  from  a  gen¬ 
eral  perspective  this  is  accurate,  the  necessary  data  processing  algorithms  are  not  de¬ 
tailed  for  the  calibration  estimations.  Hence,  suitable  estimators  need  to  be  developed 
for  each  particular  application  or  sourced  from  other  literature. 

Lassahn  and  Trenkler  (1995)  investigate  the  calibration  of  sensors  with  a  set  of  three  or¬ 
thogonal  Helmholtz  coil  pairs.  Two  methods  for  orthogonality  detection  are  presented. 
In  the  first  a  calibrated  reference  sensor  is  used  in  conjunction  with  an  interferometer 
to  detect  the  orthogonality  of  a  Helmholtz  coil  system.  Once  this  has  been  done  the 
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sensor  under  test  is  also  run  through  the  same  steps,  with  the  difference  in  field  com¬ 
ponents  indicating  the  orthogonality  errors.  In  the  second  method  a  180°  rotation  is 
used  to  cancel  out  bias  components,  providing  a  two-datapoint  measurement  for  each 
axis.  The  procedure  has  weaknesses  in  the  requirement  for  accurate  sensor  placement, 
and  physical  rotation  of  the  sensor. 

2. 1.4.1  Minimum  variance  analysis 

Auster  et  al.  (2002)  perform  a  scalar  calibration3  using  a  minimum  variance  analysis 
(MVA)  estimation.  The  authors  outline  the  parameters  which  can  be  inferred  by  rota¬ 
tions  about  a  single  axis,  or  two  or  three  axes,  proving  a  complete  parameter  estima¬ 
tion  can  be  achieved  by  performing  a  scalar  calibration  by  rotating  about  two  or  more 
axes.  The  introduction  section  also  suggests  that  it  may  be  possible  to  calibrate  a  static 
sensor,  comparing  it  to  another  stationary  calibrated  vector  magnetometer,  using  the 
variations  in  earth's  field  to  perform  the  calibration.  Whilst  this  may  be  possible,  data 
acquisition  is  most  likely  time  consuming,  and  best  done  only  for  estimating  the  rela¬ 
tionship  between  two  pre-calibrated  sensors  when  they  cannot  be  moved  due  to  their 
use,  such  as  when  using  a  sensor  as  a  reference  for  compensating  earth's  reference  field 
within  a  coil. 

2. 1.4. 2  Least  squares  fit 

Merayo  et  al.  (2000)  present  a  scalar  calibration,  using  linear  and  non-linear  least  squares 
fits.  Residual  field  strength  due  to  sensor  offsets  (bias)  is  found  to  cause  errors  for 
estimating  sensitivity.  These  errors  are  addressed  through  an  iterative  re-calibration 
leading  to  the  non-linear  least  squares  fit.  Data  needs  to  be  evenly  distributed  on  the 
surface  of  a  sphere  (representing  3D  space)  to  achieve  an  unbiased  estimate,  since  there 
is  no  weighting  to  account  for  the  angular  distribution  density  of  samples.  The  cali¬ 
bration  process  appears  resistant  to  zero  mean  Gaussian  noise.  Petrucha  and  Kaspar 
(2009)  re-implement  the  procedure  of  Merayo  et  al.  (2000),  with  three  additional  atti¬ 
tude  correction  parameters  for  the  sensor's  calibration  referencing  frame.  Additionally 

3A  scalar  calibration  uses  an  accurately  measured  total  magnetic  field  as  a  reference  signal. 
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they  extend  the  calibration  technique  to  accelerometers.  Their  results  show  the  neces¬ 
sity  of  calibrating  in  the  absence  of  any  ferromagnetic  or  high  permeability  materials, 
in  order  to  avoid  degradation  of  the  estimated  calibration  coefficients.  The  developed 
system  has  an  angular  positioning  resolution  of  only  1°,  hence  the  system  cannot  be 
used  for  vector  calibrations.  Scalar  calibrations  are  acceptable  since  the  errors  in  angu¬ 
lar  position  can  be  viewed  as  Gaussian  noise. 

2. 1.4. 3  Kalman  filtering 

Huang  and  Jing  (2008)  design  an  "Extended"  Kalman  filter,  which  characterises  the 
drag,  position  and  velocity  of  a  satellite,  and  the  orthogonality,  bias  and  sensitivity 
of  the  onboard  sensor.  The  model  generates  errors  in  the  order  of  nT  for  bias  esti¬ 
mates,  and  ~  10%  errors  in  the  sensitivity  and  orthogonality  calculations.  Similarly 
Soken  and  Hajiyev  (2011)  design  an  "Unscented"  Kalman  filter,  including  sensor  atti¬ 
tude  (but  not  sensor  orthogonality)  estimation  in  the  filter,  with  sensitivity  calculations 
separated  as  an  extension  to  the  filter.  Whilst  possibly  suitable  for  a  pico-satellite,  the 
large  variations  in  the  sensitivity  estimation  and  lack  of  orthogonality  estimation  make 
this  technique  unsuitable  for  our  purposes. 

2. 1.4. 4  TWOSTEP  maximum  likelihood  estimator 

Alonso  and  Shuster  introduce  the  TWOSTEP  algorithm  (Alonso  and  Shuster  2002a, 
Alonso  and  Shuster  2003),  initially  developed  to  estimate  bias  voltages  with  no  atti¬ 
tude  information,  but  extended  to  include  sensitivity  and  orthogonality  (Alonso  and 
Shuster  (2002c)).  To  begin  with,  a  maximum  likelihood  estimator  (MLE)  is  used  to  cal¬ 
culate  an  estimate  of  the  magnetometer  errors  using  centred  data.  This  estimate  is  then 
used  as  a  first  guess  for  a  non-centred  MLE  calculation,  allowing  the  Gauss-Newton 
method  to  select  the  global  maxima  instead  of  a  local  maxima.  The  TWOSTEP  method 
can  be  used  to  estimate  the  zero  and  first  order  terms  of  a  calibration  (bias  and  sensitiv¬ 
ity  respectively),  but  it  does  not  provide  a  computational  advantage  for  higher  order 
terms.  Fortunately,  for  the  majority  of  sensors  of  interest,  a  first  order  approximation 
of  the  response  is  all  that  is  required.  Kim  and  Bang  (2007)  use  a  genetic  algorithm 
to  generate  an  initial  estimate  of  bias  for  the  TWOSTEP  method.  Whilst  there  is  no 
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improvement  evident  in  the  final  results,  their  method  provides  an  alternative  method 
for  defining  a  first  guess  for  the  estimation  of  bias. 

Wang  (2008)  presents  a  neural  network  configuration  using  a  linearised  model  for  the 
orthogonality  correction,  assuming  fabrication  errors  of  less  than  1  degree.  The  paper 
reproduces  work  by  Olsen  et  al.  (2001)  as  an  introduction  to  scalar  calibration,  using 
a  linearised  least  squares  approach,  before  defining  the  Wang  (2008)  neural  network 
design.  Experimental  results  are  not  compared  with  the  Olsen  et  al.  (2001)  implemen¬ 
tation  or  any  other  technique;  and  with  the  final  calibrated  results  varying  up  to  55  nT, 
the  approach  seems  suitable  only  as  an  initial  estimate  of  calibration  coefficients. 

2. 1.4. 5  Ellipsoid  fitting 

Camps  et  al.  (2009)  perform  a  scalar  calibration  for  sensors,  which  minimises  the  er¬ 
rors  using  the  Levenberg-Marquadt  algorithm  to  fit  the  results  onto  an  ellipsoid.  Ini¬ 
tially  only  sensitivity  and  bias  are  modelled  (Model  1),  then  orthogonality  is  included 
(Model  2).  Results  for  the  two  models  are  compared,  performing  three  runs  for  each 
sensor,  but  are  not  compared  to  any  other  estimators.  Between  the  calibration  runs 
there  appear  to  be  a  variation  of  0.1%  in  sensitivity,  10%  in  orthogonality  and  2%  in 
bias,  but  since  the  results  are  derived  from  experimental  data  instead  of  simulation  the 
true  errors  in  the  estimation  may  be  unknown. 

Vasconcelos  et  al.  (2011)  present  the  "Geometric  Approach",  in  which  they  view  the 
sampled  data  as  the  surface  of  an  ellipsoid.  Ellipsoid  parameters  are  estimated  using 
a  MLE,  and  magnetometer  calibration  coefficients  are  then  derived  from  these  param¬ 
eters.  A  pseudo-linear  least  squares  estimate  provides  a  first  guess  of  the  parameters, 
followed  by  the  MLE  of  the  ellipsoid,  and  coefficient  extraction.  The  paper  painstak¬ 
ingly  defines  the  error  framework,  starting  with  a  description  of  the  physical  sources 
of  errors,  such  as  hard  and  soft  iron  (which  generate  permanent  and  induced  magnetic 
fields),  orthogonality  of  the  sensors,  electronic  bias  voltages,  wide-band  noise,  body 
frame  alignment  and  a  catch-all  "other"  category;  it  then  couples  the  errors  back  into 
sensitivity,  orthogonality  and  bias  errors.  Whilst  the  algorithm  does  not  provide  for 
the  alignment  of  the  sensor  due  to  the  lack  of  a  suitable  reference  frame,  a  method 
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for  the  alignment  of  the  sensor  to  a  second  vector  magnetometer  is  supplied.  Simu¬ 
lated  results  comparing  the  Geometric  Approach  to  the  TWOSTEP  method  presented 
in  Alonso  and  Shuster  (2002c)  show  good  agreement  between  the  techniques. 

Zhang  and  Gao  (2009)  present  a  model  in  which  the  sensitivity  and  orthogonality  are 
combined  into  a  single  matrix,  effectively  acting  more  as  a  cross-coupling  represen¬ 
tation.  Measurements  taken  by  performing  mechanical  sensor  rotations  in  a  uniform 
field  are  used  to  obtain  a  least-squares  estimate  of  the  parameters  for  the  ellipsoid,  de¬ 
scribed  by  the  data,  similar  to  the  estimation  method  used  by  Vasconcelos  et  al.  (2011). 
Simulated  results  show  a  variance  of  approximately  1  nT  in  calibrated  total  field,  10% 
errors  in  the  estimation  of  bias,  0.001%  error  in  the  sensitivity  and  0.5%  errors  in  the 
orthogonality  correction. 


2.1.5  Attitude  correction 

As  mentioned  above,  Vasconcelos  et  al.  (2011)  present  a  simple  attitude  correction  pro¬ 
cedure,  as  a  solution  to  the  Procrustes  problem  (Mathworks  2012).  Bracken  et  al.  (2005) 
and  Vcelak  (2006),  as  part  of  their  gradiometer  calibration,  re-align  individual  sensors 
to  the  first  sensor  in  their  gradiometer  array.  As  part  of  estimating  the  alignment  ro¬ 
tation  matrix.  Bracken  et  al.  (2005)  minimise  the  dot  product  of  each  newly  calibrated 
sensor  to  the  first  sensor.  In  the  case  of  a  scalar  calibration  this  is  obviously  not  possi¬ 
ble,  since  the  reference  sensor  only  provides  a  total  field,  but  if  it  is  assumed  that  the 
sensor  is  rotated  about  a  single  axis,  and  the  body  of  the  sensor  is  aligned  to  the  plane 
perpendicular  to  this  axis,  then  the  alignment  coefficients  can  be  obtained  to  enable 
alignment  to  the  rotation  axis.  The  estimator  would  minimise  the  variance  on  the  spin 
axis  signal  and  the  variance  of  the  sum  of  the  squares  of  the  other  two  components.  If 
additionally  the  sensor  is  rotated  on  another  orthogonal  axis,  the  complete  alignment 
matrix  to  the  sensor  body  should  be  able  to  be  derived. 
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2.1.6  Summary 

Estimation  of  the  magnetometer  calibration  coefficient  can  be  performed  by  a  range  of 
techniques,  which  generally  place  the  sensor  in  an  absolute  field  of  known  magnitude 
and  relatively  low  noise.  These  techniques  involve  rotating  the  sensor  on  various  axes, 
and  measuring  the  field  strength  detected  on  each  sensor.  Recorded  data  is  then  pro¬ 
cessed  through  a  variety  of  techniques  to  calculate  sensor  properties  such  as  sensitivity 
(nT/v),  and  orthogonality  of  sensor  outputs  to  each  other,  and  bias.  Sensors  have  been 
rotated  in  a  variety  of  ways:  with  mechanical  systems  devised  to  physically  rotate  a 
sensor;  placing  the  sensor  on  a  rotating  table;  or  using  the  motion  of  a  satellite  through 
space. 

Based  on  the  review  of  the  above-mentioned  algorithms,  the  TWOSTEP  method  by 
Alonso  and  Shuster  (2002a,  2002c)  and  the  Geometric  Approach  method  by  Vasconce- 
los  et  al.  (2011)  appear  the  most  suitable  for  implementation.  This  project  focuses  on  the 
full  implementation  of  the  TWOSTEP  and  Geometric  Approach  methods  using  MAT- 
LAB  code  in  order  to  calibrate  the  magnetic  sensors  and  magnetic  volume  generated 
within  the  Ternan  and  ELFcage  systems. 
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As  seen  in  Section  2.1.4,  magnetic  coils  are  used  to  generate  a  magnetic  field  for  cali¬ 
bration  and  excitation.  This  section  will  provide  a  brief  overview  of  a  few  coils,  their 
advantages  and  coil  calibration  techniques  found  in  literature. 

Caprari  (1995)  compares  the  uniformity  of  several  common  coil  designs:  the  Helmholtz, 
Maxwell,  Garrett,  Barker  and  Braunbek  coils.  In  each  of  these  coils  the  wire  is  wound 
into  circular  coaxial  loops,  of  varying  spacing-to-radii  ratios,  with  the  aim  of  achieving 
the  maximum  field  uniformity  possible  within  the  centre  of  the  coil  setup.  The  coil 
most  often  referred  to  in  literature  is  the  Helmholtz  coil,  in  which  the  two  circular  coils 
are  spaced  a  radius  apart.  Whilst  the  easiest  coils  to  construct  (and  the  first  invented), 
the  Helmholtz  design  also  produces  the  smallest  volume  of  homogeneity.  The  Rubens 
coil  provides  an  improvement  on  the  Helmholtz  coil,  but  should  not  be  used  where 
precise  controls  of  the  field  and  field  gradient  are  required  (Kirschvink  1992).  The 
Barker  and  Braunbek  coils,  on  the  other  hand,  produce  uniform  volumes  much  larger 
than  the  Helmholtz,  at  the  cost  of  more  complex  manufacture  and  excitation  current 
requirements. 

Square  coil  implementations  of  the  Helmholtz  design  have  also  been  built,  with  a 
slightly  larger  spacing  between  coils,  resulting  in  a  larger  volume  of  magnetic  field 
uniformity  (Kirschvink  1992). 

The  "Ternan  cos(#)"  coil  used  in  the  DSTO  MOD  system  is  based  on  the  approximation 
of  a  cos(@)  current  density  distribution  on  the  surface  of  a  cylinder.  This  concept  is 
discussed  in  detail  in  Section  3.1.1.  Specific  examples  of  this  implementation  could  not 
be  found  in  the  literature,  but  the  following  two  related  designs  were  found: 


•  Zlobin  et  al.  (2001)  present  the  design  of  a  coil  for  the  "Very  Large  Hadron  Col¬ 
lider",  based  on  Nb^Sn  superconducting  wire,  a  cos (6)  coil  and  iron  yokes.  Cross 
sectional  fields  are  very  similar  to  those  modelled  in  the  Ternan  coil.  Coils  are  fab¬ 
ricated  in  two  halves  which  are  then  epoxy-impregnated  together;  the  accuracy 
of  this  alignment  is  not  stated. 
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•  Hayes  et  al.  (1985)  presents  a  high-frequency  coil  for  use  in  magnetic  resonance 
imaging  (referred  to  as  the  "Birdcage"  resonator  coil),  in  which  the  complex 
impedance  at  a  specific  frequency  is  used  to  force  a  cos (6)  current-density  dis¬ 
tribution  on  the  coil  surface.  A  method  by  which  such  a  coil  could  be  reproduced 
for  low-frequency  use  is  also  presented  in  Section  3.1,  including  modeling  of  the 
B  field  homogeneity 

Everett  and  Osemeikhian  (1966)  describe  the  construction  of  a  "spherical  coil",  which 
is  very  closely  related  to  the  cos(0)  coil.  As  the  name  indicates,  the  "spherical  coil" 
is  wound  on  the  surface  of  a  sphere,  but  the  current  density  on  the  surface  is  still  in 
a  cos(0)  distribution.  Due  to  the  coil  being  wound  on  a  sphere  instead  of  a  cylinder, 
there  should  be  a  slight  difference  in  the  axial  field  gradient. 

The  distribution  of  the  coil  winding  in  the  cos(6>)  coils  enables  the  generation  of  a 
volume  of  magnetic  field  uniformity  that  extends  along  the  cylindrical  axis  with  height 
and  width  limited  by  the  radius  of  the  cylinder,  thereby  allowing  objects  which  had 
previously  not  been  able  to  be  inserted  into  a  magnetic  volume  to  be  tested. 


2.2.1  Coil  uniformity 

Field  uniformity  is  extremely  important  for  accurate  sensor  and  hence  coil  calibration, 
but  its  impact  on  sensors  is  often  not  well  understood.  Frix  et  al.  (1994)  compare  IEEE 
standard  calibration  coil  options  (single  coil,  circular  and  square  Helmholtz  coil  pairs), 
as  well  as  sensor  misalignment,  showing  that  the  rapid  degradation  of  the  field  unifor¬ 
mity  leads  to  several  percent  error  in  the  expected  calibration  field. 

Nissen  and  Paulsson  (1996)  investigated  the  impact  of  field  non-uniformity  of  the  coil 
on  the  calibration  of  sensors,  finding  the  deviation  from  the  true  coil  calibration  factor 
(defined  as  the  proportionality  constants  between  the  excitation  currents  and  the  in¬ 
duced  magnetic  fields),  as  a  function  of  sensor  and  coil  diameter  ratios.  Interestingly, 
for  the  Helmholtz  coil,  a  sensor  must  be  smaller  than  a  fifth  the  diameter  of  the  coils  to 
achieve  less  than  0.1%  error. 
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Carter  (1976)  associates  different  coil  design  methodologies  with  the  Taylor,  Cheby- 
shev  and  Lagrange  polynomial  approaches,  with  the  Taylor  polynomial  representing 
a  "maximally  flat"  field,  and  with  the  other  two  polynomial  methods  introducing  a 
ripple  in  the  volume  of  field  uniformity  This  approach  is  similar  to  electronic  filter 
design.  Whilst  in  many  cases  ripple  within  a  volume  may  be  an  acceptable  tradeoff 
for  larger  uniform  volume,  our  desired  magnetic  volume  may  be  used  for  gradiometer 
calibrations,  and  hence  must  have  a  maximally  flat  uniform  field. 

Bronaugh  (1995)  noted  the  lack  of  understanding  concerning  the  limitations  of  the 
Helmholtz  coil,  its  available  field  uniformity  and  the  impact  of  the  size  of  a  device  un¬ 
der  test  on  the  calibration  accuracy.  Bronaugh  needed  to  achieve  the  most  effective  use 
of  space  in  a  laboratory  environment,  similar  to  DSTO's  requirements.  Coil  calibration 
was  performed  by  Bronaugh  using  the  "ruler"  approach  (i.e.  physical  measurements 
of  coil  dimensions  are  fed  into  an  analytical  model  to  calculate  the  predicted  magnetic 
field).  A  rigorous  approach  to  verifying  other  coil  parameters  necessary  for  theoretical 
modelling,  such  as  the  number  of  turns,  is  also  provided.  Finally,  Bronaugh  also  com¬ 
ments  on  the  impact  of  loading  on  the  coil:  if  too  large  an  object  is  inserted  within  the 
coil  volume,  the  coil  excitation  current  requirements  may  be  altered,  compromising  the 
field  uniformity  due  to  mutual  inductive  coupling  to  the  object. 

2.2.2  Coil  calibration 

Drake  (1994)  presents  an  overview  of  sensor  and  coil  calibrations  within  the  National 
Physical  Laboratory  in  the  United  Kingdom.  Coil  factors  are  calculated  theoretically, 
and  validated  against  search  coils  using  a  known  20  Hz  signal.  Orthogonality  and  bias 
are  not  considered. 

Klymovych  and  Pajunpaa  (2004)  present  a  triaxial  coil  system  made  of  three  sets  of 
four  square  coils.  Coil  orthogonality  is  measured  using  a  theodolite  with  built-in 
magnetic  sensor,  presumably  placed  within  the  system  during  orthogonality  measure¬ 
ments.  Coil  factors  are  measured  by  zeroing  earth's  field  on  two  axes,  and  applying 
large  positive  and  negative  signals  on  the  third  axis.  Background  noise  is  reduced 
through  the  use  of  a  remote  reference  sensor.  Accurate  current  control  is  achieved  by 
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characterising  the  response  of  reference  resistors  to  temperature,  during  coil  excitation. 
Over  the  life  of  the  system  (years),  it  was  claimed  that  the  coil  factors  varied  by  0.01% 
and  orthogonality  varied  by  20  arc-seconds. 

Piergentili  et  al.  (2011)  design  an  earth-field  simulator  for  space  applications  using  a 
triaxial  square  Helmholtz  system.  The  system  is  validated  against  computer  simula¬ 
tions  for  the  field  uniformity  and  coil  factors,  but  orthogonality  is  not  measured,  and 
hence  not  compensated  for. 

Kuberry  (1967)  measures  coil  factors  by  zeroing  out  earth's  field  using  a  fluxgate  mag¬ 
netometer,  then  using  a  scalar  sensor  to  detect  field  strength  increases  due  to  a  subse¬ 
quent  increase  in  driving  current.  Error  bounds  arising  from  coil  and  sensor  misalign¬ 
ment,  as  well  as  sensor  bias  are  estimated  in  the  order  of  0.02%,  but  are  not  compen¬ 
sated  for.  Orthogonality  between  axes  is  not  addressed,  as  the  paper  only  examines 
the  calibration  of  single-axis  coils. 

Youguang  et  al.  (2006)  built  a  triaxial  yoke  system  for  testing  magnetic  properties  of 
materials.  Search  coils  are  used  for  field  sensing,  and  are  calibrated  following  a  manual 
alignment  in  a  solenoid  of  known  magnetic  field  strength.  Sensor-coil  misalignment  is 
considered,  and  accounted  for  by  a  rotational  transformation,  but  coil  orthogonality  is 
assumed  to  be  accurate  due  to  the  mechanical  fabrication  of  the  triaxial  yoke.  Desired 
field  strength  is  obtained  using  a  feedback  control  system.  Field  uniformity  in  the  air 
gap  does  not  appear  to  have  been  fully  considered. 

Schill  and  Hoff  (2001)  present  work  in  experimentally  measuring  the  coil  factor  for  a 
Helmholtz  coil.  The  coil  under  calibration  is  placed  within  a  larger  triaxial  coil  which 
is  used  in  conjunction  with  a  fluxgate  magnetometer  to  cancel  earth's  field  at  the  cen¬ 
tre  of  the  coil.  A  proton  resonance  magnetometer  is  then  used  to  measure  the  fields 
generated  by  the  coil  under  calibration.  This  procedure  contains  several  shortcomings 
which  Schill  and  Hoff  point  out.  Firstly,  the  earth-nulling  coil  must  be  larger  than 
the  coil  under  calibration.  Secondly,  due  to  the  temperature-dependence  of  the  wire- 
resistance  in  the  offset-coil,  changes  in  temperature  will  lead  to  a  change  in  offset  field 
strength,  and  hence  introduce  a  bias  in  the  calibration.  Finally,  the  fluxgate  sensor  used 
in  offsetting  the  earth's  field  will  introduce  a  bias  error  dependent  on  its  own  internal 
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errors.  The  temperature  dependence  of  the  coils  on  the  signal  amplitude  applied  to  the 
coils  was  measured,  and  temperatures  of  ~  25° C  to  ~  80° C  were  achieved  after  the 
coil  magnetic  field  strength  had  stabilised  for  a  set  voltage  amplitude. 

Mirzaeva  et  al.  (2012)  designed  a  single  axis  calibration  system  for  fields  of  the  order  of 
1  T.  A  step  function  is  applied  to  a  drive  coil  wound  upon  a  yoke,  and  the  increasing 
magnetic  field  is  monitored  on  a  sense  coil.  Linear  regression  analysis  is  used  to  fit  the 
coil's  response.  The  resulting  coil  factor  was  then  used  to  drive  the  system  without 
feedback,  generating  fields  within  0.76%  of  expected  values. 

Shifrin  et  al.  (2000)  present  a  four-square-coil-per-axis  system,  calibrated  with  a  scalar 
sensor.  Sensor  orthogonality  is  measured  using  methods  presented  in  the  Russian  lit¬ 
erature.  Unfortunately,  an  English  translation  of  the  paper  could  not  be  obtained  in 
order  to  provide  insight  into  the  techniques  used. 


2.2.3  Coil  limitations 

Crotti  and  Giordano  (2004)  and  Crotti  et  al.  (2006)  investigate  a  coil's  frequency  re¬ 
sponse,  and  how  it  is  impacted  by  inter-winding  inductance  and  capacitance  and  the 
skin  effect,  using  a  multi-conductor  transmission  line  approach.  Coil  factors  are  found 
to  change  by  ~  1%  between  DC  and  250  kHz,  or  less  than  0.1%  for  frequencies  below 
50  kHz.  Hence  for  our  frequencies  of  interest  (<  1  kHz)  constant  coil  factors  with  re¬ 
spect  to  frequency  can  be  assumed,  but  if  the  coils  are  driven  at  higher  frequencies,  the 
coil  factor  variation  and  the  impact  on  coil  orthogonality  should  be  addressed. 

Bronaugh  (1995)  identifies  three  primary  frequency  limiting  characteristics  of  the  coil: 


•  The  highest  frequency  limit  (and  hence  the  least  impact),  is  when  wires  are  longer 
than  ~  0.15A,  and  are  hence  comparable  to  A,  where  A  is  the  wavelength  in  the 
conductor . 

•  The  middle  frequency  limit  is  caused  by  the  induced  field  in  the  windings. 
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•  The  lowest  and  dominant  limiting  characteristic  is  the  self-resonance  of  the  coil, 
occurring  at  frequency  /o  =  2nlrLC'  w^ere  ^  and  C  are  inductance  and  capaci¬ 
tance  of  the  coil. 

Practically  the  excitation  coil  will  be  operated  well  below  resonance,  and  is  more  likely 
to  be  limited  by  the  response  of  the  driving  power  supply 
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2.3  Noise  and  Shielding 

Two  principal  noise  sources  are  likely  to  impact  the  system:  Geomagnetic  noise  and 
local  transient  noise  sources.  Geomagnetic  noise  is  spatially  correlated  over  distances 
of  several  kilometers,  due  to  the  large  physical  size  of  its  sources  (Lenz  and  Edelstein 
(2006)),  and  can  therefore  be  compensated  for  by  using  a  remote  reference  sensor  to 
feed  a  compensating  signal  to  the  coil.  Localised  transient  noise  sources  on  the  other 
hand  are  not  correlated,  and  will  require  shielding. 

2.3.1  High  magnetic  permeability  ( \i )  shields 

Mager  (1968,  1970)  presents  work  on  shielding  with  high  permeability  materials  to 
"pull"  magnetic  fields  away  from  a  specified  area,  with  greatest  effectiveness  when 
there  is  a  large  permeability  within  the  shielding  material.  Using  small  nested  shields, 
shielding  factors  of  6  x  106  have  been  measured  (Donley  et  al.  (2007)). 

Tashiro  et  al.  (2011)  experimentally  characterise  the  effect  of  shielding  materials  and 
their  impact  on  the  fields.  They  show  how  small  gaps  in  the  shield  endplates  can 
introduce  significant  fields.  The  impact  of  shielding  material  on  the  coil  factor  of  a  coil 
within  the  shield  is  also  briefly  discussed. 

Kim  et  al.  (2007)  perform  computer  modeling  on  the  effectiveness  of  y-metal  and  alu¬ 
minium  shields,  to  provide  low  reluctance  paths  for  the  magnetic  field  and  hence  gen¬ 
erating  opposing  magnetic  fields  through  induced  eddy-currents  respectively.  The  ef¬ 
fectiveness  of  the  y-metal  shield  drops  off  at  ~  10  Hz  and  above,  while  aluminium  on 
the  other  hand  begins  to  improve  in  effectiveness  for  frequencies  above  1  Hz.  A  com¬ 
bination  of  the  two  shielding  materials  allows  a  high  shielding  factor  to  be  maintained 
over  a  broad  range  of  frequencies. 

Lee  and  Romalis  (2008)  investigate  the  noise  sources  introduced  by  shielding.  Lee  and 
Romalis  (2008)  also  define  an  analytical  model  for  calculating  the  Johnson  noise  caused 
by  the  y-metal  shielding  for  a  variety  of  geometries.  Based  on  this  analytical  model, 
for  the  in-lab  shielded  system  at  DSTO,  noise  from  the  shield  amounts  to  ~  5  FT,  an 
insignificant  amount  unlikely  to  impact  on  results. 
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2.3.2  Self-shielding  coils 

Yoda  (1990)  presents  self  shielded  coils,  which  are  designed  to  reduce  the  external  field 
created  by  a  coil  system  by  placing  coils  to  counteract  its  field.  Whilst  this  could  pro¬ 
vide  reduced  coupling  to  external  structures,  there  is  no  guarantee  that  the  magnetic 
field  uniformity  within  the  system  being  shielded  is  not  compromised.  On  the  other 
hand  Boyvatn  and  Hafner  (2012)  present  a  similar  system,  using  meta-materials,  with 
significant  reduction  of  field  strength  of  ~  40  dB.  Such  a  material  would  likely  provide 
excellent  shielding  for  a  small  coil,  but  appears  limited  to  specific  frequencies. 

2.3.3  Superconducting  shields 

Aaroe  et  al.  (2009)  and  Claycomb  and  Miller  (2006)  present  superconducting  shields 
which,  whilst  providing  excellent  shielding  factors,  have  the  major  drawback  of  re¬ 
quiring  liquid  Helium  or  Nitrogen  cooling  to  attain  superconductivity. 

2.3.4  Signal  filtering 

Mengchun  et  al.  (2010)  add  digital  filtering  to  a  calibration  method  referred  to  as  a 
"double  adaptive  algorithm"  taken  from  Chinese  literature,  to  improve  performance 
in  noisy  environments.  Whilst  the  core  idea  is  commendable,  the  paper  exhibits  sev¬ 
eral  flaws  in  reasoning4,  detracting  from  the  core  idea  of  removing  noise  components 
from  the  signal  before  feeding  it  into  the  calibration  algorithm.  In  principle,  a  band¬ 
pass  filter,  fast  Fourier  transform  (FFT)  or  an  estimator  could  be  used  to  select  specific 
frequencies  which  have  low  background  noise,  using  the  alternating  signal  amplitude 
to  perform  a  sensitivity  and  orthogonality  calibration  as  a  substitute  for  the  total  field 

4Noise  is  apparently  reduced  by  applying  a  low-pass  finite  impulse  response  (FIR)  filter  to  signals 
(which  do  not  appear  to  be  adequately  pre-conditioned),  removing  much  of  the  higher  frequency  noise 
(some  of  which  was  likely  introduced  through  aliasing)  before  performing  the  calibration.  The  filtering 
provides  an  improvement  to  the  calibration,  but  the  paper  then  proceeds  to  claim  that  a  filtering  post¬ 
calibration  provides  improvements  over  pre-calibration  filtering,  presumably  mistaking  the  reduction 
in  signal  variance  due  to  filtering  for  the  reduction  in  signal  variance  caused  by  proper  calibration. 
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measurements.  This  procedure  provides  improved  signal  to  noise  ratios  (SNR),  and 
hence  a  better  estimate.  The  FFT  frequency  estimators  presented  by  Liao  (2011)  would 
serve  this  purpose.  Care  should  be  taken  when  filtering,  as  the  total  |B|  field  will  con¬ 
tain  higher  order  harmonics  of  the  AC  signal. 

2.3.5  Summary 

High  permeability  materials  used  for  magnetic  shielding  can  affect  the  field  uniformity 
within  the  volume  of  the  magnetic  field  generated  by  the  coil  system.  This  effect  is 
important  because  magnetic  shielding  is  necessary  for  sensor  calibration.  However, 
the  literature  review  failed  to  identify  any  substantial  studies  addressing  this  effect 
and  as  such,  numerical  modelling  was  conducted  to  approximately  model  the  impact 
of  high  permeability  shields  on  magnetic  field  uniformity  (Section  4). 
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3.1  Coil  Definitions  and  Derivations 

The  existing  magnetic  coil  system  from  Maribyrnong  is  modelled  as  well  as  alternate 
designs  for  comparison  of  magnetic  field  uniformity  Four  coil  designs  were  selected 
for  modelling: 


•  The  Helmholtz  coil  design  using  both  circular  and  square  coils  is  modelled  as  a 
benchmark  for  comparison  as  it  is  the  typical  coil  design  used  in  many  research 
environments. 

•  The  Barker  coil  design  is  modelled  because  a  similar  coil  design  is  in  use  at  the 
RAN  MAGFAC  facility  located  at  Steele  Point  (Sydney  Harbour),  to  quantify  the 
areas  of  field  homogeneity  ( Fleet  Command  -  RAN  Ranges  and  Assessing  Unit 
2009). 

•  The  coil  setup  at  the  Maribyrnong  facility  was  a  custom  design  for  DSTO  devel¬ 
oped  by  Dr  John  Ternan  (retired);  a  literature  search  was  conducted  to  identify 
the  exact  name  of  the  coil  design,  but  none  was  found,  as  such  the  coil  is  here 
referred  to  as  the  "Ternan  cos(F)  coil".  Modelling  is  required  to  provide  a  com¬ 
prehensive  understanding  of  the  coil's  capabilities. 

•  A  low  frequency  version  of  the  Birdcage  resonator  coil  system  (Hayes  et  ol.  1985), 
was  modelled  and  constructed  because  of  its  similarity  to  the  Ternan  cos(F)  coil 
design.  To  avoid  confusion  with  the  actual  Birdcage  resonator  design,  this  coil 
will  be  referred  to  as  the  "ELFcage"  coil  (Extremely  Low  Frequency  birdcage 
coil). 


Page  25 


3.1  Coil  Definitions  and  Derivations 


3.1.1  Coil  designs 

Coil  setups  were  modelled  in  increasing  complexity: 

•  Circular  Helmholtz  -  two  identical  circular  coils  placed  symmetrically  one  on 
each  side  of  the  experimental  area  along  a  common  axis,  and  separated  by  a  dis¬ 
tance  (2 Ls)  equal  to  the  radius  (R)  of  the  coil.  Each  coil  carries  an  equal  electrical 
current  flowing  in  the  same  direction.  See  Figure  3.1. 
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Figure  3.1:  Cross  section  of  circular  Helmholtz  coil  design,  showing  relative  spacing  of 
wires. 

•  Square  Helmholtz  -  two  identical  square  coils  (side  length  2 L)  placed  symmetri¬ 
cally  one  on  each  side  of  the  experimental  area  along  a  common  axis,  and  sepa¬ 
rated  by  a  distance  ~  0.5445  times  the  height  of  the  coil  (i.e.  Ls  =  0.5445L).  Each 
coil  carries  an  equal  electrical  current  flowing  in  the  same  direction.  See  Figure 
3.2. 

•  Barker  -  four  circular  coils  placed  symmetrically  two  on  each  side  of  the  exper¬ 
imental  area.  The  two  outer  coils  have  ~  2.26044  times  the  current  to  the  inner 
coils.  Inner  coils  are  spaced  about  ~  0.23186  times  their  radii  from  the  centre  of 
the  experimental  area,  whilst  the  outer  coils  are  ~  0.940733  times  their  radii  from 
the  centre  of  the  experimental  area  (Barker  1949,  Caprari  1995)5.  See  Figure  3.3. 

Where  Lsmau  =  0.23186R  and  Lmeif  =  0.940733R. 

5Values  are  taken  directly  from  Barker  and  Caprari  papers 
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•  Ternan  cos(6>)  -  the  Ternan  cos(6>)  coil  consists  of  Ncp  coil  pairs  spaced  uniformly 
with  respect  to  cos(F)  on  the  curved  surface  of  a  cylinder  whose  base  is  in  the 
yz  plane  and  length  is  along  the  x  axis;  each  coil  has  N;  turns,  which  is  rounded 
to  the  closest  integer.  The  same  current  is  passing  through  each  coil  turn.  At 
the  extremities  of  the  coil  the  filaments  travel  along  the  curved  boundary  of  the 
circular  base  of  the  cylinder,  in  the  yz  plane.  See  Figures  3.4  -  3.5.  The  configu¬ 
rations  of  loops  shown  in  Figures  3.4  and  3.6,  is  designed  to  provide  a  uniform 
field  within  the  cylinder  parallel  to  the  y-axis  (i.e.  By).  These  sets  of  coils  are  re¬ 
ferred  to  as  y-axis  coils.  Furthermore,  a  uniform  field  Bz,  orthogonal  to  Bt/  can  be 
generated  by  a  rotation  of  the  set  of  coil  windings  by  90  degrees  about  the  x-axis. 
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Figure  3.2:  Cross  section  of  square  Helmholtz  coil  design,  showing  relative  spacing  of 
wires. 
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Figure  3.3:  Cross  section  of  Barker  coil  design,  showing  relative  spacing  of  wires. 
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These  sets  of  coils  are  referred  to  as  z-axis  coils.  These  two  orthogonal  sets  of  coil 
pairs  therefore  generate  a  uniform  field  parallel  to  the  y  and  z  axes.  A  third  coil 
(e.g.  solenoid)  is  wound  around  the  cylinder  to  provide  a  uniform  field  along  the 
x-axis,  and  this  coil  is  referred  to  as  the  x-axis  coil. 


Figure  3.4:  Cross  section  of  Ternan  cos(F)  coil  design,  showing  relative  spacing  of  wires 
for  NCp  =  7.  Each  coil  on  the  right  hand  side  of  the  z  axis  is  associated  with  a  corre¬ 
sponding  coil  on  the  left  hand  side  of  the  z  axis,  thus  forming  a  coil  pair. 

•  ELFcage  -  the  Birdcage-resonator-like  coil  consists  of  Ncp  coil  pairs  spaced  uni¬ 
formly  with  respect  to  6  on  the  curved  surface  of  a  cylinder  whose  base  is  in  the 
yz  plane  and  length  is  along  the  x  axis,  each  coil  has  |  [Nt  sin(E)]  |  turns,  and  the 
same  current  passing  through  each  turn.  At  the  extremities  of  the  coil  the  fila¬ 
ments  travel  along  the  curved  boundary  of  the  circular  base  of  the  cylinder,  in 
the  yz  plane.  See  Figure  3.6. 

3.1.2  Approximation  of  B  for  cos(0)  coils  at  centre  of  coil 

Modelling  shows  that  the  magnetic  field  for  cos(F)  coils,  which  are  at  least  three  times 
in  length  as  in  diameter,  is  equivalent  to  the  magnetic  field  based  on  Ampere's  law  for 
an  infinitely  long  wire  given  in  equation  (3.1).  This  equation  can  be  used  to  calculate 
the  approximate  field  in  the  centre  of  the  coil  to  within  an  underestimation  of  ~  5%, 
by  summing  the  By  components  from  every  line  segment  (parallel  to  the  x-axis)  of  the 
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Figure  3.5:  Coil  filaments  as  they  travel  along  the  boundary  of  the  circular  base  of  the 
Ternan  coil. 


Figure  3.6:  Cross  section  of  ELFcage  coil  design,  showing  relative  spacing  of  wires  for 

Ncp  =  J  • 
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coil  filaments  for  each  coil  pair.  The  four  wire  filaments  that  lie  parallel  to  the  x-axis, 
and  form  a  coil  pair  are  shown  in  Figure  3.7. 


Figure  3.7:  Generic  coil  pair  for  Ternan  and  ELFcage  designs. 


Another  method  which  can  be  used  for  shorter  coils  is  obtained  by  using  the  Biot- 
Savart  law  to  calculate  the  contribution  to  the  field  from  each  loop  in  the  coil  (Cheng 
1989,  Griffiths  1999).  For  single  straight  sections  of  wire,  equation  (3.2)  can  be  used 
Cheng  (1989)  to  calculate  the  field,  a  distance  r  perpendicular  to  the  centre  of  the  wire; 
and  wire  length  2  L. 


Mo IL 


B  =  as - *7  (3.2) 

^  2nr\/L2  +  r 2 

Note  that  for  L  — >  oo  equation  (3.2)  collapses  to  (3.1).  For  the  arc  at  the  end  of  the 
cylinder,  equation  (3.3)  can  be  used  to  calculate  its  y  contribution  to  the  field  at  the 
centre  of  the  coil.  The  z  component  is  ignored  as  it  will  be  cancelled  due  to  symmetry 
Equation  (3.3)  is  derived  from  the  Biot-Savart  law  applied  to  a  section  of  a  ring  which 
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Figure  3.8:  Coil  winding  for  ELFcage  design. 


represents  the  curved  boundary  (arc)  at  the  end  of  the  cylinder.  For  each  coil  pair,  there 
are  four  arc  segments  as  shown  in  Figures  3.7  and  3.9. 

jiqILt 


B  =  a, 


2tt(L2  +  r2)3/2 


sin(0) 


(3.3) 


Figure  3.9:  Laboratory  model  of  Ternan  coil.  At  the  coil  boundary  the  wire  forms  an 
arc  on  the  surface  of  the  cylindrical  former. 

Flence  the  contribution  to  the  total  field  at  the  centre  By  from  a  single  loop  which  con¬ 
sists  of  two  arc  and  two  line  segments  is  as  shown  in  equation  (3.4),  generating  a  field 
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as  shown  in  the  right  half  of  Figure  3.7,  an  equal  contribution  is  provided  by  the  match¬ 
ing  coil  on  the  left  half  of  the  figure. 


B  =  a. 


}i0lh 


+  -)  sin(0) 


(3.4) 


Wi7+^VL2  +  r2  '  r 

For  a  series  of  Nc  single  turn  loops,  as  illustrated  in  Figure  3.8,  this  becomes  equation 


(3.5). 


B 


m )IL  ,  r  1.  . 

^WPT^{^  +  r]USmm 


(3.5) 


3. 1.2.1  Ternan  cos(0)  coil 

Since  the  Ternan  cos(0)  coil  has  the  same  number  of  turns  N/  for  each  loop,  we  multi¬ 
ply  equation  (3.5)  by  the  number  of  turns,  and  set  0 \  to  arccos(l  —  -y^py),  because  coils 
are  spaced  uniformly  with  respect  to  cos(0),  giving  equation  (3.6)  for  the  y-component 
of  the  magnetic  field  at  the  centre  of  the  coil.  The  field  would  appear  to  be  more  uni¬ 
form  if  coils  were  spaced  at  the  midpoints  of  current  positions,  but  this  introduces  the 
additional  constraint  of  Nt  being  even  to  maintain  left-right  symmetry  for  the  arc  sec¬ 
tions  of  the  loops  at  9  =  j,  or  else  the  Bz  cancelation  assumptions  leading  to  equations 
(3.3)  -  (3.5)  are  void. 


B  =  a, 


HoINtL  ,  r  1,^  .  ,  ,,  2k 

7ry/L2TT2VL2  +  r2  V  V  Nc  + 1 


)) 


(3-6) 


where  Nc  =  2 NCp,  and  NCp  is  the  number  of  coil  pairs  as  defined  in  Section  3.1.1. 


In  the  case  of  the  Ternan  cos(0)  coil  from  Maribyrnong,  an  empirical  formula  was  sup¬ 
plied,  which  takes  into  consideration  the  effects  of  the  conductors  as  they  form  a  non 
overlapping  arc  on  the  yz  plane  at  the  extremities  of  the  coil  (Clarke  2011). 


B  = 


jipKI 

r 


(3.7) 


where  K  =  22.71,  and  r  =  1.003  m 

If  we  set  I  =  1  A  in  the  case  of  NCp  =  7  and  Nt  =  6;  L  =  4  m  and  r  =  1  m: 
•  Equation  (3.7)  gives  B  =  28.4529  x  10-6  T; 
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•  Equation  (3.6)  gives  B  =  28.51651  x  10  6  T;  and 

•  MATLAB  modelling6  gives  B  =  28.51676  x  10-6  T. 

There  was  found  to  be  agreement  to  within  1%  (See  Section  7.3.1)  between  the  values 
calculated  by  modelling  (MATLAB),  and  equations  (3.6)  and  (3.7);  and  the  magnetic 
field  measured  when  the  coil  was  tested.  This  level  of  agreement  is  to  be  expected, 
since  the  winding  at  the  coil  boundary  did  not  follow  the  perfect  arc  assumed  in  (3.4), 
and  sensors  were  not  properly  calibrated.  The  above  expressions  for  the  coil  factors 
refer  to  the  y  and  z-axis  coils.  The  analogous  empirical  formula  for  the  x-axis  solenoidal 
coil  is  equation  (3.7),  with  K  =  20.6  and  r  =  0.99  m,  giving  B  =  26.148  x  10“6  T. 

3. 1.2. 2  ELFcage  coil 

Since  the  ELFcage  coil  has  the  number  of  turns  in  each  loop  related  to  cos(0),  we  mul¬ 
tiply  equation  (3.5)  by  |  [Nt  sin(0*.)]  |,  and  set  6 \  uniformly  spaced  at  intervals. 
|  [Nt  sin(0/f)]  |,  is  rounded  to  the  nearest  integer,  representing  the  number  of  turns  at  0/c. 


B  =  a. 


yoiL 


Nc 


nV L2  +  r2  L2  +  r2 


+  -)L  sin( 


irk 


k=l 


Nc  + 1 


[Nt  sin( 


nk 


Nc  + 1 


)]| 


(3.8) 


where  Nc  =  2 NCp,  and  Ncp  are  defined  in  Section  3.1.1. 


*cp 


For  large  Ncp  the  Ternan  cos(0)  coil  approaches  a  perfect  cos (6)  current  distribution; 
and  if  Nt  is  sufficiently  large  so  will  the  ELFcage,  as  such  it  can  be  stated  that  the  Ternan 
and  ELFcage  coils  are  simply  two  different  approximations  of  the  cos($)  coil. 


3.1.3  Engineering  considerations 

Selection  of  Ncp  and  Nt  values  will  affect  cost,  ease  of  fabrication,  as  well  as  the  uni¬ 
formity  of  the  magnetic  field  within  the  coil.  Intuitively  the  most  uniform  field  will  be 
achieved  using  a  Ternan  cos($)  coil  with  an  large  value  of  Ncp  and  Nt  =  1,  as  it  will 
provide  the  closest  approximation  to  the  field  generated  by  a  cylinder  with  a  cos(0) 
current  distribution  along  its  surface. 

6Full  coil  was  modelled  using  Biot-Savart  law  to  compute  magnetic  field  throughout  the  volume 
within  the  coil. 
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3.1  Coil  Definitions  and  Derivations 


Figure  3.10:  Laboratory  model  of  Ternan  coil. 


A  major  disadvantage  to  this  approach  is  that  the  construction  of  the  coils  becomes 
quite  complicated  and  prone  to  errors  in  the  placement  of  wires,  due  to  strict  require¬ 
ments  in  their  positioning,  and  the  sheer  number  of  separate  loops  to  be  wound  (See 
Figures  3.9,  3.10);  hence  wires  need  to  be  grouped  together  either  by  increasing  Nt 
whilst  decreasing  Ncp,  or  by  moving  to  the  ELFcage  design  with  a  lower  value  of  Ncp, 
and  higher  Nt-  This  complication  is  compounded  by  the  addition  of  a  complete  set  of 
orthogonal  coil  pairs  required  to  produce  a  magnetic  field  uniformity  parallel  to  the 
z-axis. 

The  the  estimated  wire  requirements  outlined  in  equations  (3.9)  and  (3.10).  The  Ternan 
cos(F)  coil  in  the  DSTO  magnetic  test  facility  has  values:  Ncp  =  7  and  Nt  =  6;  L  =  4 
m  and  r  =  1  m.  The  ELFcage  was  designed  with  Ncp  =  7  and  Nt  =  14;  L  =  4  m  and 
r  =  1  m.  Hence  the  Ternan  and  ELFcage  coils  require  ~  lAkm  and  ~  2.1  km  of  wire 
respectively  The  differing  wire  requirements  do  not  directly  impact  on  manufacture 
cost,  since  differing  wire  gauges  would  be  required,  but  the  process  of  winding  an 
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ELFcage  would  most  likely  be  more  time-consuming  than  the  winding  of  a  Ternan 
coil. 


I 


W Ternan 


Nc/ 2  2k 

2 nr  +  4 LNcNt  +  8 Ntr  ^  arccos(l  -  — - -) 

k=i  Nc  +  1 


(3-9) 
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W  ELF cage 
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fc=i  iVc  ~l~ 1  fc=i 


nk 

Nc  +  1 


)] 


nk 

Nc  +  1V 


(3.10) 


The  advantage  of  the  Ternan  cos(0)  coil  with  a  lower  NCp,  but  higher  N;  is  that  fabrica¬ 
tion  becomes  much  easier:  all  loops  have  the  same  number  of  turns,  are  spaced  further 
apart  than  for  large  NCp,  and  the  maximum  field  strength  increases  linearly  with  N); 
but  uniformity  is  lost  near  the  boundaries  of  the  coil  along  the  y-axis. 

The  advantage  of  the  ELFcage  coil  that  for  the  same  value  of  NCp  and  sufficiewntly 
large  Nt,  a  better  uniformity  is  achieved  near  the  y-axis,  due  to  the  uniform  spacing  of 
coils  with  respect  to  9.  Uniformity  very  close  to  the  boundary  is  reduced  compared  to 
a  perfect  cos(0)  coil,  but  not  as  much  as  for  the  Ternan  approximation. 
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3.2  Simulation 

3.2.1  Volume  calculations 

MATLAB  code  was  written  to  reproduce  each  of  the  modelled  coils  as  a  series  of 
straight  line  segments  defining  the  path  the  wire  would  follow  when  actually  wind¬ 
ing  the  coil.  In  the  case  of  the  Helmholtz,  Barker  and  Ternan  cos(6>)  coils,  the  in¬ 
creased  number  of  turns  in  the  coil  was  modelled  by  a  linear  increase  of  current  passing 
through  a  single  wire;  whereas  in  the  case  of  the  ELFcage  coil  each  winding  was  indi¬ 
vidually  included.  The  BiotSavart  law  was  then  used  to  calculate  the  contribution  from 
each  segment  to  our  area  of  interest. 

3. 2. 1.1  Field  strength  at  centre  of  coil 

From  the  field  strengths  Bx,  By  and  Bz  calculated  in  MATFAB,  Bcentre  was  defined  to  be 
the  field  strength  at  the  centre  of  the  coil  as  shown  in  equation  (3.11).  Bx  and  Bz  were 
found  to  be  zero  at  the  centre  of  the  coil. 


B 


centre  — 


B 


}) centre 


(3.11) 


3. 2. 1.2  Field  normalisation 

Since  coils  were  designed  to  generate  a  field  in  the  y  direction,  fields  were  normalised 
as  shown  in  Table  3.1,  to  provide  a  clear  representation  of  the  variation  of  the  field  with 
respect  to  the  centre. 

Results  for  the  cross  section  at  the  centre  of  each  coil  are  in  Figures  3.11,  3.12, 3.13, 3.14 
and  3.15,  where  all  variations  greater  than  5%  have  been  clipped  to  allow  visibility  of 
small  changes. 
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Field  Component  axis  Normalisation 


Reason 


A  B 


X% 


100 1 
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B centre  ' 


Calculate  %  variation  from  zero 


A  B 


10Q||By|  iWel- 


y% 


A  B 


2% 


100|J^-|fc 

'  Drpntrp  ' 


Calculate  %  variation  from  B, 


centre 


Calculate  %  variation  from  zero 


Table  3.1:  Field  normalisations,  where  i,  j,  k  are  unit  vectors  along  the  x,  y,  z  axes. 
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Note  that  in  Figures  3.11  -  3.15,  in  the  lower  right  panel,  at  the  furthest  points  of  the  coil 
in  the  x-z  plane,  the  coils  appear  connected  by  a  wire  parallel  to  the  y  axis.  This  wire 
is  in  fact  a  pair  of  wires  with  current  running  in  opposite  directions,  and  was  included 
to  faithfully  represent  how  the  coils  are  wound  in  real  life.  The  wire-pair  does  not 
contribute  to  the  total  magnetic  field  as  the  wires  cancel  out  each  other's  field. 

3. 2. 1.3  Volume  surface 

Data  points  were  considered  to  be  part  of  the  uniform  volume  if  they  met  three  condi¬ 
tions: 


•  ABtotaix  was  less  than  desired  uniformity  error  bounds  (see  equation  (3.12)). 
Equation  (3.12)  represents  the  absolute  value  of  the  sum  of  the  column  vectors 
ABX%,  A By%,  A Bz%. 

•  The  point  was  within  a  radius  0.85r  from  the  x-axis,  to  reject  the  filaments  of 
ABtotai%  ~  0  field  error  which  appear  between  closely  spaced  wires.  These  fila¬ 
ments  are  clearly  shown  as  the  ripple  in  the  field  uniformity  in  Figures  3.14  and 
3.15. 

•  Starting  from  the  point's  projection  on  the  x-axis,  all  points  to  the  centre  of  the 
coil  along  the  x-axis  are  within  the  volume. 


ABtotalo/o  —  |  &BX%  +  ABy%  +  A  Bz%\  (3.12) 

Figures  3.16  -  3.20  show  the  ±1%  A Bt0taio/  volumes  for  the  Helmholtz,  Barker,  Ternan 
cos(F)  and  ELFcage  coils;  Figures  3.21  -  3.25  show  the  ±3%  equivalents.  What  appear 
to  be  steps  along  the  x-axis,  are  in  fact  caused  by  the  discrete  number  of  samples  taken 
along  said  axis. 
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Y  Axis  (m)  -1  -1  X  Axis  (m) 


Figure  3.16:  Volume  bounded  by  a  ±1%  variation  from  Bcentre  within  Square  Helmholtz 
coil. 


Figure  3.17:  Volume  bounded  by  a  ±1%  variation  from  Bcentre  within  Helmholtz  coil. 
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Figure  3.18:  Volume  bounded  by  a  ±1%  variation  from  Bcentre  within  Barker  coil. 
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Figure  3.19:  Volume  bounded  by  a  ±1%  variation  from  Bcentre  within  Ternan  cos(F)  coil. 
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Figure  3.20:  Volume  bounded  by  a  ±1%  variation  from  Bcentre  within  ELFcage  coil. 


3.2  Simulation 


Figure  3.21:  Volume  bounded  by  a  ±3%  variation  from  B centre  within  Square  Helmholtz 
coil. 


Figure  3.22:  Volume  bounded  by  a  ±3%  variation  from  Bcentre  within  Helmholtz  coil. 
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Figure  3.23:  Volume  bounded  by  a  ±3%  variation  from  Bcentre  within  Barker  coil. 
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Figure  3.25:  Volume  bounded  by  a  ±3%  variation  from  Bcentre  within  ELFcage  coil. 


3.2  Simulation 


3.2.2  Comparison  of  coil  designs 

As  a  measure  for  comparing  the  effectiveness  of  each  coil,  the  ratios  of  the  uniform 
volume  to  the  total  volume  bounded  by  the  coil  and  of  the  uniform  area  to  the  total 
area  in  the  yz-plane  are  shown  in  Figures  3.26  and  3.27  respectively.  Usable  volume 
and  area  were  calculated  using  the  method  defined  in  Section  3.2. 1.3. 


Usable  uniform  volume 


Figure  3.26:  Comparison  of  usable  volume  as  a  function  of  the  percentage  field  unifor¬ 
mity  (ABtotal%  (3.12):  0  -  10%). 


As  shown  in  Figure  3.26,  the  Ternan  cos(O)  and  ELFcage  coils  provide  a  large  uniform 
volume  to  work  with  for  any  set  uniformity  requirements,  with  reasonably  steep  in¬ 
creases  in  usable  volume  as  the  requirements  are  relaxed.  The  ELFcage  exhibits  strong 
radial  growth  at  the  centre  as  values  are  relaxed  from  0  to  1%  (See  Figure  3.27),  making 
it  ideal  for  ranging  of  short  objects  with  large  diameters  to  a  high  degree  of  accuracy. 
The  ELFcage  and  Ternan  cos(F)  coil  are  superior  to  the  Helmholtz  and  Barker  designs 
only  under  specific  uniformity  values.  Their  performance  could  be  improved  by  an 
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increase  of  NCp  for  both  coils,  and  an  increase  in  Nt  for  the  ELFcage.  Preliminary  sim¬ 
ulations  indicate  the  ELFcage  can  generate  a  more  uniform  field  than  the  Barker  coil 
when  both  NCp  and  Nt  are  increased  to  values  of  ~  10  and  ~  16  respectively. 

Note  that  the  Barker  coil  appears  to  be  the  most  suitable  given  that  it  has  the  highest 
usable  volume  and  area  ratios.  Flowever  the  Barker  coil  has  a  lower  diameter  to  length 
ratio  compared  to  the  cos(0)  coils  for  creating  a  uniform  field.  For  example  for  3%  field 
uniformity  Figure  3.23  shows  that  2  m  diameter  coils  are  required  to  generate  a  uniform 
By  field  that  extends  approximately  1  m  in  length  along  the  x-axis.  Therefore  a  longer 
prolate  object,  e.g.  6  m  in  length,  would  require  a  coil  diameter  of  approximately  12 
m,  which  is  impractical.  However,  Figures  3.24  and  3.25  show  that  the  field  uniformity 
for  the  Ternan  and  ELFcage  coils  extends  to  approximately  6  m. 


Usable  uniform  area  on  YZ  plane,  at  centre  of  coil 


Figure  3.27:  Comparison  of  usable  area  as  a  function  of  the  percentage  field  uniformity 
(A Btotai%  (3.12):  0  -  10%). 
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The  "Finite  Element  Method  Magnetics  4.2"  package  (Meeker  2010),  is  used  to  investi¬ 
gate  the  impact  /i-metal  shielding  would  have  on  the  field  uniformity  of  a  small  triaxial 
magnetic  coil  system. 


A  solenoidal  field  was  modeled  in  free  air,  inside  a  cylindrical  shield  with  one  open 
end,  and  inside  a  closed  cylindrical  shield  with  a  small  gap  in  the  junction  between  the 
body  of  the  shield  and  the  lid.  An  infinitely  long  ELFcage  coil  was  also  modeled  with 
and  without  shielding. 


4.1  Solenoid 


The  solenoid  has  the  following  physical  dimensions:  radius  0.28  m,  length  0.66  m, 
with  340  turns  of  2  mm  diameter  wire,  having  0.1  A  flowing  through  them.  The  fi- 
metal  shield  is  in  the  form  of  a  can,  and  is  2  mm  thick  with  an  internal  radius  of  0.320 
m,  and  is  0.684  m  deep.  The  lid  provides  an  additional  0.058  m  of  space  above  the 
shield  body.  The  lid  is  2  mm  thick,  and  is  modeled  as  having  an  outside  diameter  that 
is  6  mm  larger  than  the  outside  diameter  of  the  shield  body,  leading  to  a  1  mm  gap  in 
the  overhang  between  the  lid  and  the  body  (See  Figure  4.1),  this  gap  is  a  worst-case 
scenario  for  incorrect  lid  placement,  and  in  real  life  is  expected  to  be  of  the  order  of  0.1 
mm. 
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646mm 


684mm 


Not  to  Scale 


Figure  4.1:  Illustration  of  fz-metal  shielding  container. 


Figures  4.2,  4.3,  and  4.4  illustrate  how  the  magnetic  field  lines  are  pulled  toward  the 
shielding  material,  with  angles  of  incidence  close  to  90°,  and  the  magnetic  field  lines  in 
Figure  4.4  are  visibly  straightened  by  the  end-fire  shielding  from  the  lid.  Note  that  the 
contour  lines  are  not  an  indication  of  field  strength,  since  they  were  plotted  using  an 
auto-scale  function.  Figures  4.3,  4.5  and  4.6  show  the  effect  of  removing  the  shielding 
lid  (located  at  approximately  —0.4  m  in  Figures  4.5  and  4.6),  causing  the  field  at  the 
coil  boundary  to  be  pulled  radially  outward  toward  the  shield  body,  creating  a  greater 
field-gradient  than  in  the  in-air  coil.  The  cross  sections  in  Figures  4.7,  and  4.8  show  the 
extent  to  which  field  uniformity  is  broadened  by  the  gradual  increase  in  shielding. 

Figure  4.5,  shows  the  effect  of  the  gap  between  the  shield  lid  and  the  coil  (at  approx¬ 
imately  —  0.4  m)  compared  to  the  shield  base  and  the  coil  (at  approximately  +0.4  m), 
leading  to  a  slight  reduction  in  field  strength.  Compare  right  and  left  hand  sides  of  the 
"green"  curve  in  Figure  4.5. 

The  spike  in  magnetic  field  strength  at  approximately  +0.4  m  (Figure  4.5)  from  the  cen¬ 
tre  is  due  to  the  boundary  between  air  and  the  shield,  followed  by  a  the  final  extracted 
data-point  being  set  to  zero  for  plotting  purposes;  this  glitch  is  also  present  in  Figure 
4.7,  and  to  a  lesser  extent  in  the  other  figures.  The  asymmetry  in  the  in  air  coil  in  Figure 
4.5  ("red"  curve)  is  caused  by  the  coil  and  shield  not  being  co-centred,  that  is  the  field 
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lines  are  computed  to  the  edge  of  the  can  on  the  right  hand  side,  and  to  the  edge  of 
the  lid  on  the  left  hand  side,  hence  causing  the  air-gap  to  the  left  of  the  solenoid  to  be 
included. 

Shielding  will  tend  to  increase  the  uniformity  of  the  generated  field. 


Figure  4.2:  Magnetic  field  generated  by  solenoid  in  air.  Diagram  shows  right  hand  side 
of  field  lines.  The  leftmost  line  represents  the  solenoid  axis  of  symmetry. 


Figure  4.3:  Magnetic  field  generated  by  solenoid  in  open  ^ -metal  container.  Diagram 
shows  right  hand  side  of  field  lines.  The  leftmost  line  represents  the  solenoid  axis  of 
symmetry. 
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Figure  4.4:  Magnetic  field  generated  by  solenoid  in  closed  y-metal  container.  Diagram 
shows  right  hand  side  of  field  lines.  The  leftmost  line  represents  the  solenoid  axis  of 
symmetry. 

4.1.1  Anomalies 

Figures  4.6  and  4.8  exhibit  a  sharp  drop  in  the  field  at  the  origin,  this  is  due  to  FEMM 
not  calculating  the  field  at  the  simulation  boundary,  whence  when  interpolating  the 
cross  section  plot  near  the  boundary  a  straight  line  is  generated  linking  the  closest 
data-point  within  the  finite  element  mesh  and  the  boundary  As  the  mesh  density  is 
increased7,  the  drop  becomes  steeper,  since  the  data-points  are  closer  to  the  origin. 


7 


The  mesh  resolution  was  increased  as  more  details  were  added  to  the  geometry 
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x  iq-5  Solenoid  |B|  Field  Strengths 


Figure  4.5:  Total  field  strength  along  solenoid  axis  of  rotational  symmetry.  The  lid  is 
located  at  —0.4  m,  and  the  base  is  located  at  +0.4  m. 


x  iq-5  Solenoid  |B|  Field  Strengths 


Figure  4.6:  Total  field  strength  along  the  radial  direction  from  the  axial  centre  of 
solenoid. 
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Normalised  Solenoid  |B|  field 


Figure  4.7:  Total  field  strength  variation  along  solenoid  axis  of  rotational  symmetry. 
The  lid  is  located  at  —0.4  m,  and  the  base  is  located  at  +0.4  m. 


Normalised  Solenoid  |B|  field 


Figure  4.8:  Total  field  strength  variation  along  the  radial  direction  from  the  axial  centre 
of  solenoid. 


Page  60 


Chapter  4 


Effect  of  Magnetic  Shielding  on  Field  Uniformity 


4.2  ELFcage 


The  simulated  ELFcage  had  the  following  physical  properties:  Shield  radius  is  0.32  m, 
ELFcage  radius  is  0.28  m,  Nj  =  11,  and  NCp  =  4,  with  a  current  of  0.01  A  running 
through  the  coil  filaments. 

Figures  4.9  and  4.10  seem  to  indicate  the  shielding  does  not  significantly  improve  uni¬ 
formity  within  the  ELFcage  coil,  however  it  does  show  qualitatively  the  increase  in 
field  strength  with  shielding  compared  to  an  in-air  coil.  In  Figure  4.9  the  blue  circu¬ 
lar  boundary  represent  the  fictitious  boundary  of  the  shield.  In  Figure  4.10  the  blue 
circular  boundary  represent  the  boundary  of  the  shield. 

Similar  to  the  case  of  the  shielded  solenoid.  Figures  4.11  and  4.12  show  how  the  mag¬ 
netic  field  strength  is  amplified  by  the  presence  of  shielding;  the  generated  field  strength 
is  increased  by  a  factor  of  two,  while  the  induced  currents  in  the  shields  slightly  im¬ 
prove  uniformity,  as  shown  in  Figures  4.13  and  4.14.  Figures  4.11  to  4.14  represent  the 
fields  along  the  two  orthogonal  axes  matching  the  y  and  z  axes  in  Figure  3.7  in  Section 
3.1.2. 

The  currents  in  the  shield  are  only  a  slightly  better  approximation  to  a  cos  6  current 
density  distribution  than  the  wires,  hence  if  a  better  approximation  of  a  cos  6  current 
density  coil  (large  NCp )  were  simulated,  the  shield  would  only  negligibly  improve  the 
magnetic  field  uniformity. 
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Figure  4.9:  Magnetic  field  generated  by  infinitely  long  ELFcage  in  air. 


Figure  4.10:  Magnetic  field  generated  by  infinitely  long  shielded  ELFcage. 

4.3  Summary 


Modelling  shows  a  pi- metal  shield  will  increase  both  the  uniformity  and  the  field  strength 
of  solenoidal  and  ELFcage  coils.  In  solenoidal  coils  the  relative  distance  between  the 
circular  shield  lids/endcaps  to  the  coil's  end-fire  boundary  compared  to  the  shield 
body  to  coil  end-fire  boundary  has  a  strong  impact  on  field  uniformity  For  best  unifor¬ 
mity,  the  shield  end-caps  must  be  as  close  as  possible  to  the  solenoid's  end-fire  bound¬ 
ary 

A  finite  length  ELFcage  has  field-distorting  currents  running  at  the  coil  boundaries,  in 
the  same  location  as  and  flowing  in  the  same  plane  as  the  solenoid  end-fire  boundary 
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1Q-''  ELFcage  |B|  Field  Strengths 


Figure  4.11:  Total  field  strength  along  horizontal  (y)  axis  of  ELFcage. 


1Q-''  ELFcage  ]B|  Field  Strengths 


Figure  4.12:  Total  field  strength  along  vertical  (z)  axis  of  ELFcage. 
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Normalised  ELFcage  |B|  field 


Figure  4.13:  Normalised  total  field  strength  variation  along  horizontal  (y)  axis  of 
ELFcage. 


Normalised  ELFcage  |B|  field 


Figure  4.14:  Normalised  total  field  strength  variation  along  vertical  (z)  axis  of  ELFcage. 
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Shield  end-plates  are  expected  to  pull  associated  field-lines  away  from  the  main  body 
of  the  ELFcage,  significantly  improving  field  uniformity 

High  permeability  shielding  is  expected  to  improve  magnetic  field  uniformity  for  solenoid 
and  ELFcage  coil  designs,  whilst  concurrently  attenuating  external  magnetic  noise, 
hence  //-metal  shielding  is  ideal  for  isolating  laboratory  setups  in  noisy  environments. 
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5.1  TWOSTEP 


This  section  "walks  through"  the  derivation  of  the  complete  TWOSTEP  method  in 
greater  detail  than  presented  in  Alonso  and  Shuster  (2002c).  Latter  sections  address  re¬ 
defining  the  coil  and  sensor-within-coil  calibration  models  so  as  to  match  the  TWOSTEP 
algorithm,  enabling  it  to  be  used  to  estimate  their  properties.  For  ease  of  comparison 
to  the  Alonso  and  Shuster  papers,  cross-references  to  the  relevant  papers  are  included. 
References  [{x}]  and  [(x)]  indicate  the  Xth  equation  as  marked  in  Alonso  and  Shuster 
(2002a)  and  Alonso  and  Shuster  (2002c)  respectively,  where  both  papers  contain  the 
same  equation  both  will  be  referred  to,  as  [{x},(x)]. 


5.1.1  Model 

We  define  the  sensor  model  as  follows  [(45)]: 

Bk=(I  +  D)-'(#TAtHk+b  +  £t)  (5.1) 

where  B %  is  the  kt]l  sensor  reading,  H \  is  earth's  field  at  that  reading,  b  is  sensor  bias,  e ^ 
is  the  measurement  noise.  D  is  a  fully  populated  symmetric  matrix  representing  sen¬ 
sitivity  and  orthogonality  errors  when  added  to  the  identity  matrix  I.  A/(  is  an  attitude 
correction  matrix  (rotation)  used  to  put  H/c  in  the  same  reference  frame  as  B/(,  and  i9  is  a 
rotational  error  within  the  sensor,  hence  another  rotation  matrix.  Since  only  |ff^|  is  ex¬ 
pected  to  be  known  during  the  calibration  (A^  and  d  are  not  known  or  estimated),  (5.1) 
is  manipulated  to  eliminate  the  rotations  (which  are  singular  orthogonal  matrices),  and 
use  only  magnitude  data  [(19)]. 
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Thus: 


(I +  D)Bk  =  dT  AkHk  +  b  +  ek  (5.2a) 
(7  +  D)Bk  -  b  -  ek  =  VTAkHk  (5.2b) 

((7  +  D)Bk  -  b  -  ek)T((I  +  D)Bk  -  b  -  ek)  =  \Hk\2  (5.2c) 

5.1.2  Parameter  re-definition 

I  Hk  | 2  is  re-arranged  and  a  better  set  of  basis  functions  is  selected  for  improved  model 
stability  (Barker  et  al.  2007),  hence: 

I Hk\2  =  ((7  +  D)Bk  -  b  -  ek)T((I  +  D)Bk  -  b  -  ek)  (5.3a) 

=  Bk  (7  +  D)t(7  +  D)Bk  -  Bk  (7  +  D)Tb  -  BTk  (7  +  D)Tek 
-bT(I  +  D)Bk  +  bTb  +  bT£k  (5.3b) 

—  £k  (I  +  D) Bk  +  £kb  +  £k£k 
=  Bl (7  +  2D  +  D2)Bk  -  2 BTkc  -  2 BJk  (7  +  D)T£k 

(5.3c) 

+  2bT£/f  +  \£k\2  +  |b|2 
=  |B,.|2  +  BlTEBt-2Bfc 

(5.3d) 

-(2[(J+D)Bt-b]-£t-|£t|2)  +  |b|2 

=  | B,  | 2  +  hIfRi  -  2Bjc  -  vt  +  |b|2  (5.3e) 

Where  for  ease  of  estimation  we  removed  the  non-linear  dependence  on  parameter  D, 
defining  matrix  E  and  vector  c  as  shown  in  (5.4)  [(46)],  and  the  scalar  vk  is  explicitly 
defined  in  (5.7b). 

E  =  2D  +  D2  (5.4a) 

c  =  (7  +  D)b  (5.4b) 

As  we  can  see  c  is  simply  the  bias  recast  through  sensitivity  and  orthogonality  errors. 

Manipulating  (5.4)  as  shown  in  (5.5)  and  (5.6)  expresses  b  in  terms  of  c  and  E,  therefore 
b  =  b(c,  E),  (5.6c). 
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2D  +  D2 

(5.5a) 

c  =  (J+D)b 

(5.6a) 

2J  +  D 

(5.5b) 

(i  +  oy 

_1c  =  b 

(5.6b) 

I  +  D 

(5.5c) 

1 

h-1 

tn 

1 

xc  =  b 

(5.6c) 

5.1.3  Sensor  “error  measurement”  &  “error  measurement  noise” 

We  work  with  scalar  measurement  zk  and  scalar  measurement  noise  vk,  as  defined 
in  [(23a,  24)];  which  is  the  nomenclature  used  in  the  Alonso  and  Shuster  papers.  zk 
and  vk  should  not  be  confused  with  ek  and  Bk  being  measurement  noise  and  the  /cth 
sensor  measurement  respectively.  As  shown  in  (5.3d)  and  (5.3e),  v k  holds  all  the  e k 
noise  elements  present  in  \Hk\2,  and  zk  could  be  considered  a  scalar  measurement  of 
the  sensor  errors  we  are  trying  to  estimate,  zk  and  vk  are  defined  in  (5.7). 


zk  =  \Bk\2-\Hk\2  (5.7a) 

vk  =  2[(I  +  D)Bk-b]-ek-\ek\2  (5.7b) 


Using  (5.3e)  the  scalar  measurement  zk  is  expressed  as  follows  [(47)]: 

(5.8a) 
(5.8b) 


zk  =  \Bk\2  -  \Hk\2 

=  -Bj.EBk  +  2Blc  +  vk-  \b\2 
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5.1.4  Vectorisation 


For  ease  of  derivation  and  coding,  E  and  c  are  vectorised,  rearranging  to  remove  the 
symmetric  3x3  matrix  E  and  to  use  6x1  vector  E  instead,  as  done  in  (5.9)  [(48)]. 


/  \  /  r„  ,  \ 


En  E12  Ei3 
E 12  E  22  E  23 


Elk 
Eik 

\  B3k  J 


Bk  EBk  -  y  Blk  Blk  B3k  , 

\  E\3  E2 3  E33 
—  EnB>\k  +  E2iB2k  +  E33B3k 
+  Ei2(2BlkB2k)  +  Ei3(2  BlkB3k)  +  E23{2B2kB3k) 


—  {  B\k  B2fc  B*k  2BlkB2k  2BlkB3k  2B2kB3k 


=  KiE 


(  En  ^ 
En 
£33 
En 

E13 

V E23 ) 


(5.9a) 


(5.9b) 


(5.9c) 


(5.9d) 


Where  Kk  [(49a)]  and  E  [(49b)]  are  defined  as  follows: 


E  —  En  E22  E33  E12  Ei3  E23 
Kk  =  (  B2fc  B2fc  B2^  2BlkB2k  2B3kB3k  2 B2kB3k 


(5.10a) 

(5.10b) 


Applying  this  simplification  to  zk  in  (5.11)  [(50)]  and  merging  E  and  c,  zk  is  expressed 
in  terms  of  variables  Lk  and  d'  shown  in  (5.12)  [(51)]. 


—  BkEBk  +  2Bkc  +  vk  —  b|3 

(5.11a) 

—KkE  +  2Bkc  +  vk—  |b(c,  E)|2 
/  \ 

(5.11b) 

(  2B\  -Kk  )  (  C  )  —  | b ( c,  E )  | 2  +  vk 

(5.11c) 

v  v E  / 

Lt9'~  |b(0')|2  +  i/t 

(5. lid) 
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L/c  is  1  x  9  row  vector  dependent  on  sensor  measurements  Bk/  and  0'  is  a  9  x  1  column 
vector  dependent  on  the  redefined  sensor  errors  E  and  c,  hence  on  D  and  b. 


Lk  = 


(5.12a) 


0'  = 


(5.12b) 


The  use  of  0'  as  opposed  to  0,  reflects  the  need  to  recover  b  and  D  from  0'.  0  = 
(  b  D  )T  can  be  derived  from  0'  as  is  shown  later  in  Section  5.1.10.  A  re-arrangement 
of  (5.11)  represents  vk  in  terms  of  Lk  and  0'  as  shown  in  (5.13). 

vk  =  Zk-Lk0'  +\b(0')\2  (5.13) 


5.1.5  Noise  distribution 

As  part  of  the  MLE  process  the  noise  in  the  system  needs  to  be  characterised.  ek  is 
assumed  to  be  white  and  Gaussian  distributed  with  zero  mean.  Since  ek  r-u  M(o,zk) 
[{5},(4)]  then  14  ~  as  shown  in  (5.14a)  [{7a},(5a)]  and  (5.15a)  [{7b},(5b)]. 


Calculation  of  }ik  is  shown  below: 

n  =  E  M 

(5.14a) 

=  £{2[(J+D)Bt-b]  £t}-E{|et|2} 

(5.14b) 

=  0-e{(^a^)2} 

(5.14c) 

=  -fr(Et) 

(5.14d) 

f  4,i  o  o  ) 

0  el 


where  E*.  is  the  covariance  matrix  for  sample  k  is 


0 


,  and  fr(Efc) 
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Calculation  of  cr2  requires  the  use  of  Gaussian  random  variable  moments  presented  in 


Appendix  A.l. 

=  E  {vk}  -  }4  (5.15a) 

=  4E  {([(/  +  D)Bk  -  b]T£t)2}  -  2E  {([(/  +  D)Bk  -  b]T£t)|£t|2}  +  E  {l£l|4}  -  tr  ( Lk  f 

(5.15b) 

=  4E  {[(/  +  D)Bk  -  b]T£t£tT[(l  +  D)Bk  -  b]}  +  E  {|£l|4}  -  tr  (Zk)2  (5.15c) 

=  4[( I  +  D)Bk  -  b]TZt[( I  +  D)Bk  -  b]  +  51  r  (z2)  -  3(r  (z2)  (5.15d) 

=  4[( I  +  D)Bk  - b]TZt[( I  +  D)Bk-b]+  21  r  (zf  )  (5.15e) 

where  3b(L2)  =  [fr(E^)]2  (See  equation  (A.13)). 


5.1.6  Likelihood  function 

Having  characterised  the  noise,  the  likelihood  (5.16a)  and  negative  log-likelihood  (5.17a) 
[(6)]  functions  are  obtained: 


N 

/n|0,=6/(0/) = n 

fc= i 


1  p-(vk-Vk)2 

i 


N 


=n 


1  c-(zk-Lke' +\b(e')\2-uky  /2a^ 


(5.16a) 

(5.16b) 


/  (&')  -  —^nfN\0'=0'(O') 


-1  N 
z  k= 1 


1  2 
^2  (Z^  -  M'  +  |b(e')  I2  -  +  In  (27rcrf ) 


(5.17a) 

(5.17b) 

(5.17c) 
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5.1.7  Centering 

As  part  of  the  TWOSTEP  method,  J  (df)  will  be  split  into  a  centred  and  centre  compo¬ 
nent,  keeping  the  non-linear  components  in  the  centre  component.  (5.18),  and  (5.19) 
define  the  centre  variables  and  (5.20)  defines  the  centred  variables  [(7,  8,  9,  52)]. 


1  N  1 

u  k= 1  ak 

N| 

III 

1-1 

N 

?r* 

N  - 1  N  ■ | 

B  =  i r2E^Bk  L  =  a2  Y  -jL, 

k=  1  ak  k= 1  iTk 

N 

v  =  v1  E  -2v* 

k=  1  ak 

N  ^ 

k= 1  ak 

Zk  =  Zk-Z 

Bk  =  Bk  —  B  Lk  =  Lk  —  L 

Vk  =  Vk-  V 

/A  =  ft 

Lk6'  +  vk  and  z  =  Ld' 

-  b(0')  |2  +  v. 

(5.18) 


(5.19) 


(5.20) 


/  ( 6 ')  can  be  split  into  centre  and  centred  terms  as  shown  in  (5.21)  [(10)]. 


ii 

© 

1 

(5.21a) 

1  N  -t 

5  E[^(v/c-^)2  +  ln27ro-f] 
z  k= 1  Pfc 

(5.21b) 

1  N  i 

5  E  [^(^  +  V  -  Hie  -  H)2  +  ln  27nrfc2] 

Z  fc=l  ak 

(5.21c) 

1  N  i 

oEb((h-  h)2  +  iy-  v)1  +  2(vjfc  -  fik)(v  -  p))  +  ln27rc^] 

z  fc=i  ^ 

(5.21d) 

i  N  i 

/  (d')  +  /  (0')  +  x  E  —  flk)(v  ~  p)] 

Z  fc=l 

(5.21e) 

J  (8')  +  J  (8') 

(5.21f) 
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Where  /  (5.22)  is  the  centred  measurement  component  and  /  (5.23)  is  the  centre  mea- 

N  ^ 

surement  component  [(56,  58)],  and  ^  —2(2(14  —  flk)(v  ~  fi))  =  0  as  shown  in  Ap¬ 


pendix  A. 2. 


k=l  ak 


1 


N 


J  {O')  =  3  E  -  /ifc)2  +  ln2 nal 

z  fc=i 

1  N  1  ~  2 

=  -  7 1  —  (zfc  —  Lfrd'  —  fa)  +  terms  independent  of 

2  fc=l 


(5.22a) 

(5.22b) 


1 


N 


/(0/)  = 

Z  /c=l 

=  2^2  (2-^0/+|b(0')|2-H' 


(5.23a) 

(5.23b) 


5. 1.7.1  Centered  estimator 


Estimator  for  the  centred  measurement  is  found  by  setting  = 0  [(57a)]- 

0  0  ( \  N  \ 

J  =  ——  -  ~2  (zfc  —  Lkd'  —  fib)  +  terms  independent  of  6'  J  (5.24a) 


W  d0'\2  k-ak 


1  N  —2Lt 

E  — 2^  (2*  ~  ^0'  -  n) 


2fc“  ak 


N 


-1 


H  lit 


Z-^(2k-n)LTk  +  E^?0'  =  ° 


N  fT 


UL 


k=  1  °"k 

N  -1 


fc=l 


E  -t10'  =  E  ^2  (Z*  -  fo)  U 


k=l  ak 


k=  1  ^ 


(5.24b) 

(5.24c) 

(5.24d) 


N  fT 


LtL 


We  multiply  both  sides  by  the  inverse  of  ^  -  2  /v  to  obtain  the  centred  estimate  0' 


fc=l 


[(57a)],  where  *  indicates  the  value  is  an  estimate. 


1 


,  N  TtT,\  1  N 

r=(Ei^  L^_{it_h)L 

4=1  ak  J  k= 1  ak 


N  1 


fc=i  ^ 


(5.25a) 

(5.25b) 


Where  5  (d/)  is  the  centred  Fisher  Information  Matrix  and  is  derived  in  Section  5.1.8. 
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5.1.8  Fisher  information  matrix 


The  Fisher  Information  Matrix  3  is  the  expectation  of  the  Hessian  matrix  -  the  square 
matrix  of  the  second  order  partial  derivatives  of  the  negative  log  likelihood  function 
/  (6').  The  fisher  matrix  is  defined  by  Poor  (1994)  as  (5.26): 


3  (0')  —  Ex|©'=0'{( 


dln/X|©/=e/(X)  91n/X|0/=0,(X) 


dd' 


-)(- 


dd' 


(5.26) 


The  Fisher  matrix  can  also  be  represented  as  (5.27),  as  shown  in  White  (2012). 


3  (»')  =  -Ex|e-=0-{(3  ln/3y2#(X)»  (5.27) 

5. 1.8.1  Centered  Fisher  Information  Matrix 

Using  (5.27)  the  Fisher  Information  Matrix  for  /  (0')  is  derived8,  as  shown  in  (5.28) 
[(57b)]. 


3  (o')  =  EU^ljpP-n  (528a) 

Tp  i  N  i 

=  E{(  w  9  E  -  Lk6'  -  FkV(zk  -  Lk0'  -  \ik ))}  (5.28b) 

00  z  fc=i  a k 

d  N  -l 

=  E{^L^2LTk(^-L^'-fik)}  (5-28c) 

°°  fc=i  ak 
N  -1 

=  £{E  -Uh}  (5.28d) 

/c— 1  ak 
N  i 

=  £  -^ll~Lk  (5.28e) 

A:=l  ak 

5. 1.8. 2  Center  Fisher  Information  Matrix 

As  stated  in  (5.26)  and  (5.27),  there  are  two  ways  of  calculating  the  Fisher  matrix,  for 
the  centre  log-likelihood  function  /  (0')  [(59)],  we  will  use  (5.26). 

8Due  to  the  differentiation  5  (i/ j  is  no  longer  dependent  on  9'.  9'  is  nonetheless  maintained  in  the 
notation  to  distinguish  3  (‘;)  from  3  (‘) 
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3(0') 


EX|0=fl{( 


a/(0').T,a/(0') 


30' 


)} 


1 


Ex|©=e{(^2)  + 


30' 

a|b(9')|2 

30' 


)T(y-fl)' r(v-p)(-I  + 


,M|2 


5|b(0')l 

30' 


02(-I  + 


i/\|2 


3|b(0')| 


V2'  v  30' 

i/r  3|b(e')lVrr 

^2l  30'  M 


Ex|0=e{(v  _  + 


/M2 


5|b(9')l 

30' 


A|2 


9|b(0')l 

30' 


(5.29a) 

)} 

(5.29b) 

(5.29c) 

(5.29d) 


Where  EX|@=0{(v  —  ^)9}  is  derived  in  Appendix  A.3,  and  the  derivative  of  |b(0')|2 
with  respect  to  0'  is  shown  below. 

|b(0')  |2  is  presented  in  terms  of  c  and  E  [(54)]: 


|b(0')|2  =  b(0')Tb(0') 

=  cT((I  +  D)~1)t(I  +  D)~1c 
=  ct(I  +  D)~2c 
=  ct(7  +  2D  +  D2)“1c 
=  cT(I  +  E)~1c 


(5.30a) 

(5.30b) 

(5.30c) 

(5.30d) 

(5.30e) 


Splitting  0'  into  sub-vectors: 


_3_ 

S07 


|b(0')|2 


l\He')\2 

^|b(0')|2 


(5.31a) 


|b(0')  | 2  is  then  differentiated  element  by  element  with  respect  to  c  and  E  [(55)]: 


A|b(0')|^Acr(I  +  E)-lc 

(5.32a) 

=  ((I+E)-1c  +  ct(I+E)-1), 

(5.32b) 

=  (2(J  +  £)-1c)l 

(5.32c) 
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m 

(.hi) 

E  element 

E  element 

O'  element 

1 

(1/1) 

En 

Ei 

04 

2 

(2,2) 

E22 

Ei 

05 

3 

(3,3) 

E33 

e3 

06 

4 

(1,2) 

E11 

e4 

07 

5 

(1,3) 

E13 

e5 

08 

6 

(2,3) 

E13 

Ee 

09 

Table  5.1:  Mapping  E  indexed  elements  'Em'  and  matrix  E  elements 


3EJb(fl')'2  =  3E,/T<'  +  E)-1c 

(5.33a) 

(5.33b) 

~  ct(/  +  E)-'3('  +  ^(/  +  E)-'c 

(5.33c) 

W|b(0')|2  =  -(2-<S,v)(cT(/  +  E)-1)i((J  +  £)-1c);  (5.34a) 

=  -(2  -  >,/)((!  +  E)~1c)i((I  +  E)~1c)j  (5.34b) 

Where  5^  is  the  Kronecker  delta  function,  and  i  &  j  take  on  the  values  shown  in  Table 
5.1. 

Consideration  of  (5.21)  and  (5.27)  implies  the  complete  Fisher  matrix  can  be  obtained 
by  adding  the  centre  and  centred  components  as  in  (5.35). 

3  (0')  =  3  (0')  +  3  (0')  (5.35) 
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5.1.9  Center  term  correction 

The  Gauss-Newton  method  is  used  to  iteratively  bring  the  estimate  O'*  closer  to  the 
true  value  O'.  Starting  with  O'q  =  O'* ,  we  iterate  (5.36)  until  ;/,  in  (5.37)  reaches  a 
predetermined  small  quantity,  or  i  reaches  an  upper  limit  [(14a,  14b)]. 

0'+1=8'-3(0')-1d7/(0')  (5.36) 

m  =  W  -  0'-i)T3(0'-l)  W  -  S'-i)  (5.37) 


5. 1.9.1  Log-likelihood  gradient 

Section  5.1.9  required  the  gradient  of  /  (O')  to  be  evaluated  at  O'*.  Equation  (5.21)  is 
used  to  derive  the  ^  gradient  by  evaluating  components  ^  and  ^  separately  ^  / 
was  derived  in  (5.24)  and  refined  in  (5.38).  is  fully  derived  in  (5.39). 


Pi  N  _  I  N  TT  f 

^I=L4  (4-  -  M  LI  +  E  (5.38a) 

ou  k=l  Uk  k=  1  ak 

N  i 

=  E  —  (4  -  F‘k)  LTk  +  5  (e')e'  (5.38b) 

k=  1  ak 

r)  3  1  N  1 

^J=Mkl2L-f-P)2]  (559a) 

=  +  lb(0,)f  “  P)  1  (5.39b) 

=  -1  (ie'-  ^7lb(0')|2)  (z-i0'+|b(0')|2-)i)  (5.39c) 
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5.1.10  Converting  to  b  and  D 

Estimate  O'*  is  in  terms  of  E*  and  c* ,  as  opposed  to  the  desired  D*  and  b*.  E*  can  be 
changed  back  into  matrix  form  E*  by  using  its  definition  in  (5.10a).  To  compute  D*  we 
write  [(60)]: 


E*  =  USUT  (5.40) 

where  U  is  orthogonal  and  S  is  a  3  x  3  diagonal  matrix  with  elements  s,  [(61)],  and  both 
matrices  are  obtained  by  the  Eigen-decomposition  of  symmetric  matrix  E*. 

We  define  diagonal  matrix  W,  satisfying  S  =  2 W  +  W2  [(62)],  hence  using  the  quadratic 
formula  the  elements  are9  [(63)]: 

Wj  =  —  1  T-  yj  1  +  s j  (5.41) 

The  maximum  likelihood  estimate  for  D  can  hence  be  calculated  by  [(64)]: 

D*  =  UWUT  (5.42) 

Bias  vector  b*  can  be  calculated  as  shown  in  (5.6c)  [(65)]: 

b*  =  (I  +  D*)-V  (5.43) 

5.1.10.1  Covariance  Matrix 

To  obtain  the  covariance  matrix  we  need  to  transform  3  (9')  to  3  (0),  where  6  = 

and  D  is  defined  with  respect  to  D,  the  same  way  as  E  was  defined  in  (5.10a)  [(67)].  Ex¬ 
plicit  derivation  of  this  process  is  not  shown  [(66,  68,  69)]. 


3(e)'1  = 


(5.44a) 

(5.44b) 


9Wj  —  —  1  —  y/1  —  Sj,  is  not  used  as  it  would  give  negative  sensitivity  values,  implying  the  sensor  was 
facing  the  wrong  way.  This  is  not  possible  since  it  would  have  been  accounted  for  by  attitude  rotations. 


and  manufacture  quality  control  is  expected  to  capture  such  major  sensitivity  errors 
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Element  by  element  differentiation  of 


3(b,D) 

9(c,E) 


can  be  represented  as 


fw)  = 

V  rW  )  i,j  dei 


Mb  ,D)\  = 

(/+D)  McD(b) 
06x3  21  +  Med(D) 


(5.45a) 

(5.45b) 


with 


and 


M£D(D) 


b\ 

0 

0  b2 

b3 
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McD(  b)  = 

0 

b2 

0  bi 

0 

b3 
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(  2D1 
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2D  2 
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0 
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1  +  C>2 

d6 
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d5 

0 

d5 

D6 

D 

+  C>3  D4 

V  0 

d6 

d6 

d5 

d4 

d2  +  d3 

\ 


(5.46) 


(5.47) 


5.1.11  Implementation 


Begin  by  assuming  |H/C|,  B/c,  and  are  known. 


•  Calculate  p^  =  —  fr(Ejt) 

•  Calculate  z ^  using  (5.7a). 

•  Calculate  L ^  using  (5.12) 

•  Set  6'  =  0,  and  calculate  i using  (5.15a). 

•  Calculate  centre  components  a1,  B,  z,  and  p  using  (5.18),  and  (5.19). 

•  Calculate  centred  components  B ^  z k,  L^,  and  pi,  using  (5.20). 

•  Calculate  centred  Fisher  Information  Matrix  ^(f^*)  using  (5.28). 
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•  Estimate  O'*  using  (5.25). 

•  Calculate  log  likelihood  gradients  using  (5.38)  and  (5.39). 

•  Calculate  centre  Fisher  Information  Matrix  5(0'*)  using  (5.29). 

•  Iterate  centre  term  correction  as  described  in  Section  5.1.9,  updating  5(0'*)  at  each 
iteration. 

•  Recover  desired  values  b  and  D  as  shown  in  Sections  5.1.10  and  5.1.10.1. 

An  optimised  implementation  of  the  code  can  be  found  in  Appendix  C.l.  I  gratefully 
acknowledge  Dr  Jose  Vasconcelos  for  providing  his  implementation  of  the  TWOSTEP 
code  for  debugging  and  validating  this  implementation. 
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5.2  Geometric  Approach  -  Ellipsoid  Calibration 


This  section  briefly  summarises  the  model  and  estimator  presented  in  Vasconcelos  et  al. 
(2011),  along  with  a  minor  correction  to  their  first-guess  estimation.  The  model  was 
implemented  for  comparison  with  TWOSTEP,  but  not  taken  beyond  the  calibration  of 
sensors  in  constant  magnetic  fields,  as  the  model  struggled  with  large  data-sets,  even 
when  code  was  optimised. 


5.2.1  Model 

The  sensor  model  is  defined  as  follows: 


hri  =  CBhi  +  b  +  nmi  (5.48) 

where  /z„  is  the  zth  vector  reading  from  the  sensor;  C  =  SmCnoCsi  II  Eh  ||  is  a 
fully-populated  matrix  built  from  diagonal  scaling  matrix  Sjvf/  lower-triangular  non¬ 
orthogonality  matrix  Cjvo,  and  fully  populated  soft-iron  transformation  matrix  Cg  / ; 
b  =  SmCno^hi  +  is  sensor  bias,  with  hard-iron  bias  bni  and  null-shift  1?m;  nmi  is 
noise  on  the  zth  reading;  and  Bhj  =  BREhi  is  the  true  magnetic  field  rotated  from  the 
earth's  field  reference  frame  to  the  sensor-body  frame. 

Vasconcelos  et  al.  (2011)  show  the  calibration  parameters  can  be  inferred  from  the  cen¬ 
tre,  orientation  and  radii  a  3D  ellipsoid,  the  surface  of  which  is  defined  by  the  best  fit 
of  acquired  sensor  data. 

Initially  a  least  squares  estimation  is  used  to  obtain  a  first  guess  p,  a  9  x  1  vector,  which 
is  the  solution  to  the  least  squares  problem  (5.49): 


1 

* 
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h2 
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where  p  = 


i  T 


A  B  C  D  E  G  H  I  J 


and  hrix,  hriy,  and  hnz  are  the  zth  x,  y,  and 


z  components  of  recorded  data  hn.  p  is  converted  to  the  initial  estimate  of  calibration 
parameters  Tq10  (5.50)  and  bo  (5.51): 


Tn  = 


~l  tan(p) 

■i(tan(p)  tan(A)  sec (cp)  —  tan (cp)) 


+  lsec(p) 


0 

0 


—  I  sec  (p)  tan(A)  sec(cp)  ~c  sec(A)  sec  (cp) 

(5.50) 


bo=2^ 


h 

p2 

h 


(5.51) 


Parameters  a,  b,  c,  p,  A  and  cp  are  numerically  calculated  in  Appendix  A.4.1,  and  defined 
in  Section  5.3. 


In  the  second  stage  the  vectorised* 11  initial  estimate  xq  = 


bo 


calculated  above  is 


used  to  place  the  first  guess  of  the  MLE  close  to  the  estimator's  peak.  The  Gauss- 
Newton  approach  is  used  to  bring  the  solution  close  to  the  true  maxima: 


-  (V2/Ut)-‘V/Ut  (5.52) 

where  V/ \Xk  and  V2/|.Y;c  are  the  gradient  and  Hessian  matrices  of  the  log-likelihood 
function  f(x),  and  are  re-produced  in  A.4.2. 


Tk+l 

Tk 

bk+i 

- 1 

A 

10  The  sign  of  the  term  at  the  centre  of  the  Tq  matrix  has  been  changed  from  that  presented  in  Vascon- 
celos  et  al.  (2011). 

11 21  is  the  vector  defined  by  the  stacked  column  elements  of  T/,. 
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5.3  Estimator  Selection 

In  this  section  processing  time  and  accuracy  of  the  Geometric  Approach  and  TWOSTEP 
estimators  are  compared,  and  one  is  selected  for  further  use.  The  key  measures  of 
effectiveness  will  be  accuracy,  processing  time,  estimator  stability  and  algorithm  per¬ 
formance  for  large  data-sets.  Un-calibrated  "measured"  and  "true"  data  is  initially 
compared  with  their  calibrated  total  field  measurement  for  data-sets  of  80  and  8000 
samples.  Performance  of  the  two  estimators  then  is  tested  over  a  range  of  noise  levels 
and  sample  sizes.  Throughout  all  the  simulations  the  sensor  is  stimulated  from  ran¬ 
dom  directions,  describing  the  surface  of  a  sphere,  as  is  illustrated  in  Figure  5.1.  A 
Geometric  Approach  calibration  on  pre-calibrated  TWOSTEP  data  is  also  performed, 
in  an  attempt  to  reduce  total  processing  time  when  using  Geometric  Approach.  This 
method  is  referred  to  the  "TS  +  GA"  method. 

Stimulation  data,  as  illustrated  in  Figure  5.1  was  generated  using  the  MATLAB  "rand" 
function  to  assign  random  8  and  tp  spherical  coordinates  to  the  magnetic  field  direction, 
for  a  constant  magnitude  of  one,  forming  the  surface  of  a  sphere.  Normally  distributed 
noise  is  then  added  to  the  magnetic  field  amplitude  of  each  axis,  with  noise  variance 
dependent  on  simulation  parameters,  consequently  moving  data-points  above  and  be¬ 
low  the  surface  of  the  sphere. 

Referring  to  Figure  5.1,  note  that: 

•  Assuming  perfect  orthogonality,  good  sensitivity  equates  to  the  magnitude  of  the 
semi-principal  axes  of  the  measurement  ellipsoid  being  equal  to  each  other  and 
also  to  the  magnitude  of  the  field  component  being  measured,  thereby  creating  a 
unit  sphere,  which  may  be  displaced  from  the  origin,  dependent  on  bias. 

•  Assuming  perfect  sensitivity,  good  sensor  orthogonality  equates  to  alignment  of 
the  measurement  ellipsoid's  axes  to  the  cartesian  axes.  A  non-orthogonality  of 
the  sensor  will  cause  the  measurement  ellipsoid  axes  to  be  rotated,  and  a  slight 
deformation  from  the  shape  of  a  sphere. 

•  Good  bias  equates  to  the  origin  of  the  measurement  ellipsoid  and  the  sphere  be¬ 
ing  coincident. 
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Calibrated  and  Uncalibrated  data 

*  T^true 

*  B measured 


Figure  5.1:  Simulated  "true"  (green)  and  "measured"  (blue)  magnetic  field  components 
of  B. 


Large  sensitivity  (Kx,  Ky  and  Kz)  and  orthogonality  (p,  A  and  (p)  errors  were  selected  to 
represent  a  badly  built  sensor  or  coil,  and  bias  (bX/  by  and  bz)  was  included  in  the  simu¬ 
lation,  but  is  not  presented  in  the  results  since  the  Geometric  Approach  and  TWOSTEP 
models  place  it  in  differing  locations,  of  their  respective  models,  meaning  that  moving 
either  bias  from  one  model's  framework  to  the  other's  would  couple  sensitivity  and 
orthogonality  estimation  errors  into  the  resulting  bias,  thereby  influencing  the  result; 
hence  the  RMS  value  of  |B  +  Noise \  —  \R  Estimated]  is  used  instead.  Numerical  values  of 
the  errors  are  shown  in  Table  5.2. 
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Kx 

Ky 

Kz 

P 

A 

<P 

bx 

by 

bz 

3.9589 

2.9243 

2.2485 

-2.0608 

-19.4292 

16.5846 

0.0055 

0.0137 

0.0235 

Table  5.2:  Measurement  errors  used  in  the  simulation.  Units  of  Kx,  KXJ,  Kz,  bx,  by,  and 
bz  are  arbitrary,  and  p,  A,  cp  are  in  units  of  degrees. 

5.3.1  First  order  performance 

The  "true"  total  field  strength  in  the  following  section  is  unity  with  isotropic  noise  of 
variance  of  1.25  x  1CU3  added  to  the  dataset12.  Total  fields  for  two  datasets  of  80  and 
8000  samples  are  compared  to  get  a  first  order  understanding  of  the  performance  of 
the  estimators. 

Figures  5.2  and  5.3  show  the  true,  un-calibrated  and  calibrated  total  field  |B  +  noise \ 
data,  for  the  TWOSTEP  and  Geometric  Approach  estimates.  The  data  calibrated  using 
the  two  estimates  closely  follows  the  simulated  total  field  and  magnetic  noise,  along 
with  what  appears  to  be  a  slight  estimator  bias  in  the  sensitivity  and  bias  error  esti¬ 
mates,  leading  to  an  underestimate  of  the  total  field.  As  sample  size  is  increased  the 
total-field  bias  appears  to  reduce,  as  would  be  expected  since  MLEs  are  asymptotically 
unbiased. 

Figures  5.4  and  5.5  show  the  difference  between  |B  +  noise |  and  |Besf;-mate|.  The  estima¬ 
tor  bias  alluded  to  above  is  evident.  The  TWOSTEP  bias  appears  to  be  slightly  larger 
than  that  in  the  Geometric  Approach,  with  both  approaching  zero  as  sample  size  is 
increased. 

5.3.2  Reference  frame  selection 

TWOSTEP  and  Geometric  Approach  impose  different  restrictions  on  the  calibration 
matrices,  which  impacts  on  the  way  sensitivity  orthogonality  errors  are  represented. 
The  sensitivity  and  orthogonality  correction  in  the  TWOSTEP  process  can  be  expressed 
symbolically  as  shown  in  equation  (5.53). 

H  =  KB  (5.53) 

12£i  —  £2  —  £3  —  1-25  x  1  (3  3 .  Measurements  are  relative  to  each  other,  not  absolute. 
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Calibrated  and  Uncalibrated  data 


Figure  5.2:  Comparison  between  uncalibrated  data  with  noise  (pale  blue)  and  true 
(known  magnetic  field  data)  with  noise  (green),  calibrated  (estimated)  estimate]  using 
TWOSTEP  (blue)  and  Geometric  Approach  (red)  methods.  80  Samples. 


T  T 

where  H  =  [Hx,  HXJ,  Hz]  is  the  calibrated  signal,  B  =  [Bx,  By,  Bz]  is  the  raw  mea- 

Kr  a  b 


sured  signal,  and  K  = 


a  Ky  c 
b  c  Kz 


is  the  calibration  matrix. 


Hj  =  KXBX  T  ciBy  T  bBz  (5.54a) 

Hy  =  aBx  +  KXjBy  +  cBz  (5.54b) 

Hz  =  bBx  +  cBy  +  Kz  Bz  (5.54c) 


TWOSTEP  uses  a  fully  symmetric  matrix,  which  effectively  means  each  component 
from  the  raw  data  e.g.  Bx  contributes  an  equal  proportion  of  itself,  to  each  of  the  other 
two  calibrated  axes  Hy  and  Hz(via  the  off-diagonal  elements),  as  the  corresponding 
raw  data  for  each  of  said  calibrated  axes  contributes  as  a  proportion  of  themselves  to 
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Figure  5.3:  Comparison  between  uncalibrated  data  with  noise  (pale  blue)  and  true 
(known  magnetic  field  data)  with  noise  (green),  calibrated  (estimated)  \Bestimate\  using 
TWOSTEP  (blue)  and  Geometric  Approach  (red)  methods.  8000  Samples. 
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Figure  5.4:  TWOSTEP  and  Geometric  Approach  |B  +  noise  |  Errors  -  80  Samples. 
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Figure  5.5:  TWOSTEP  and  Geometric  Approach  |B  +  noise |  Errors  -  8000  Samples. 


the  original  calibrated  axis.  Equation  (5.54)  illustrates  this,  where  Hx  receives  a  from 
By  and  b  from  BZr  and  Bx  contributes  a  to  Hy  and  b  to  Hz. 

The  Geometric  Approach  uses  a  lower  triangular  matrix  as  a  first  guess,  meaning  the 
x-axis  is  assumed  to  have  no  orthogonality  error,  and  the  error  projected  onto  the  y  and 
z  axes;  but  the  final  Geometric  Approach  matrix  has  no  restrictions  on  it,  meaning  there 
are  no  restrictions  on  the  reference  frame  of  the  sensor.  To  enable  the  comparison  of 
specific  calibration  parameters,  results  must  be  brought  into  a  single  reference  frame. 

For  convenience  a  lower  triangular  matrix,  was  selected  as  the  reference  frame  (shown 
in  Figure  5.6)  for  comparing  results,  this  reference  frame  matches  frame  used  in  the 
first  guess  in  the  Geometric  Approach.  TWOSTEP  and  Geometric  Approach  sensitiv¬ 
ity  and  orthogonality  matrices  are  easily  brought  into  the  correct  frame  using  QL  de¬ 
composition,  which  produces  a  singular  orthogonal  matrix  Q  (a  rotation  matrix)  and  a 
calibration  matrix  L  shown  in  Equation  (5.55),  from  which  sensitivity  ( Kx ,  Ky  and  Kz) 
and  orthogonality  (p,  A  and  (p)  errors  can  be  extracted  algebraically.  L  can  be  used  to 
transform  the  ellipsoid  (x  ,  y  ,  z  axes.  Figure  5.6)  in  Figure  5.1  to  the  sphere  (x,  y,  z  axes. 
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Figure  5.6)  representing  the  true  field,  by  multiplying  the  ellipsoid  data  by  L,  using  the 
values  of  Kx,  Kx/,  Kz,  p,  A,  and  <p  from  Table  5.2. 


Figure  5.6:  Coordinate  framework  for  lower-triangular  matrix,  where  (  x 
T 


V  2 


L 


v 


converts  the  raw  measured  sensor  readings  to  calibrated  readings. 


Kx  0  0 

KXJsin(p)  KXJcos(p )  0 

Kz(sin((p)  cos(A))  Kzsin(A.)  Kzcos(cp)  cos(A) 


(5.55) 


5.3.3  Performance  for  noisy  data 

The  total  field  strength  in  the  following  section  is  unity  with  a  isotropic  noise  of  vari¬ 
ance  of  2  x  ICC6  to  2  x  ICC2  added  to  the  data-set  of  2000  samples.  Performance  char¬ 
acteristics  were  measured  at  each  data-set  and  averaged  over  30  simulations  with  the 
same  noise  characteristics.  Note  the  spikes  in  the  Geometric  Approach  curves  occur 
when  the  calibration  algorithm  has  failed  to  converge  to  a  solution  within  acceptable 
stop  conditions1^.  Use  of  the  median  of  the  simulation  runs  would  have  removed  the 
spikes  and  provided  a  better  approximation  to  the  expected  performance,  but  it  was 

13The  stop  conditions  relate  to  the  magnitude  of  the  step  size  at  each  iteration. 
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believed  important  to  retain  the  anomalies  as  an  illustration  of  the  impact  certain  data¬ 
sets  can  have  on  the  estimator.  Estimator  stopping  conditions  were  set  to  the  point 
where  solutions  just  start  not  being  able  to  converge,  to  enable  a  fair  comparison  of  the 
best  achievable  solutions  for  each  algorithm. 


5.3.3. 1  Processing  Time  and  Number  of  Iterations 

Figure  5.7  shows  the  processing  time  for  the  Geometric  Approach  is  slightly  depen¬ 
dent  on  noise  variance,  with  a  slight  increase  in  high  noise  environments,  whilst  the 
TWOSTEP  algorithm  exhibits  a  minima  at  10-4  noise  variance,  rapidly  increasing  if 
noise  variance  is  reduced,  or  slowly  increasing  if  noise  is  increased.  The  reason  the 
TWOSTEP  processing  time  increases  for  low  noise  is  associated  with  the  numerical 
precision  in  MATLAB,  leading  to  quasi  singular  or  badly  scales  matrices,  when  a  small 
noise  variance  is  used.  This  is  most  likely  caused  by  the  centring  process  (Section  5.1.7). 
Estimation  using  TWOSTEP  is  significantly  faster  than  using  Geometric  Approach. 


Processing  Time  vs  Noise  Variance 


Figure  5.7:  TWOSTEP  (blue).  Geometric  Approach  (red),  and  TS  +  GA  (magenta): 
processing  time  vs  noise. 
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The  reason  for  the  Geometric  Approach's  close-to-constant  processing  time  becomes 
apparent  from  Figure  5.8.  The  Geometric  Approach  tends  to  converge  in  approxi¬ 
mately  5  iterations,  unless  the  data  is  pre-processed  with  TWOSTEP  ,  in  which  case 
it  converges  in  as  little  as  3  iterations.  TWOSTEP  on  the  other  hand  takes  between  5 
and  100  iterations  to  converge  to  a  solution. 


Figure  5.8:  Number  of  TWOSTEP  (blue).  Geometric  Approach  (red),  and  TS  +  GA 
(magenta):  iterations  required  for  convergence  vs  noise. 


5. 3. 3. 2  \B  +  noise\  RMS  error 

Figure  5.9  shows  the  RMS  error  of  the  total  field  for  both  estimators  as  noise  is  in¬ 
creased.  TWOSTEP  and  Geometric  Approach  are  effectively  identical  in  their  perfor¬ 
mance  for  low  noise,  but  Geometric  Approach  marginally  outperforms  TWOSTEP  for 
high  noise  environments. 
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Figure  5.9:  TWOSTEP  (blue).  Geometric  Approach  (red),  and  TS  +  GA  (magenta): 
RMS  of  (|B  +  noise \  —  estimate])  vs  noise. 


5. 3. 3. 3  Individual  parameter  errors 

For  a  direct  comparison  between  TWOSTEP  and  Geometric  Approach  calibration  pa¬ 
rameters,  both  sensitivity  and  orthogonality  estimates  were  extracted  using  QT  decom¬ 
position,  as  shown  in  Section  5.3.2. 

Figure  5.10  shows  that  the  TWOSTEP  and  Geometric  Approach  orthogonality  errors 
increase  with  increasing  noise.  Geometric  Approach  appears  to  slightly  outperform 
TWOSTEP  for  high  noise  environments,  but  fails  to  converge  more  often  than  TWOSTEP, 
causing  in  "spikes"  in  the  simulation  results. 

Similarly  the  sensitivity  errors  shown  in  5.11  show  good  agreement  between  the  two 
methods,  but  the  Geometric  Approach  clearly  gets  slightly  better  in  high  noise. 
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Figure  5.10:  TWOSTEP  and  Geometric  Approach:  angle  errors  vs  noise. 


Figure  5.11:  TWOSTEP  and  Geometric  Approach:  sensitivity  errors  vs  noise. 
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5.3.4  Performance  for  large  data-sets 

The  total  field  strength  in  the  following  section  is  unity  with  an  isotropic  noise  of  vari¬ 
ance  of  10“4  added  to  the  data-sets  consisting  of  20  to  2  x  105  samples.  Performance 
characteristics  (processing  time  and  number  of  iterations)  and  estimation  errors  were 
recorded  for  each  data-set  and  averaged  over  30  iterations  with  the  same  noise  charac¬ 
teristics. 

5.3.4. 1  Processing  time  and  Number  of  Iterations 

Figure  5.12  shows  the  increase  in  processing  time  as  sample  size  is  increased.  TWOSTEP 
outperforms  Geometric  Approach  by  between  one  and  two  orders  of  magnitude.  When 
optimising  the  code  it  was  noticed  the  bulk  of  the  time  within  Geometric  Approach  is 
taken  up  by  Kronecker  product  (<g>)  operations  over  k  samples  and  to  a  lesser  extent 
normal  matrix  multiplication  operations  between  matrices  that  are  k  samples  deep14. 
This  slow-down  is  caused  by  the  large  memory  requirements  necessary  to  process  all 
the  data  concurrently  in  the  matrices,  and  is  likely  due  to  the  calibration  matrix  not 
being  triangular  or  symmetric,  hence  not  being  vectorised  to  reduce  memory  usage. 
Extrapolating  from  the  graph,  107  samples  (2hours  45min  of  data  sampled  at  1kHz) 
could  be  processed  by  TWOSTEP  in  10s. 

Figure  5.13  indicates  processing  time  is  principally  dependent  on  memory  usage,  since 
the  number  of  iterations  for  the  Geometric  Approach  is  fundamentally  flat.  TWOSTEP 
on  the  other  hand  seems  to  dramatically  increase  the  number  of  iterations  it  takes  to 
converge  once  the  sample  size  increases  above  2000  samples. 

5. 3. 4. 2  \B  +  noise\  RMS  error 

Figure  5.14  shows  the  RMS  error  for  the  total  field  is  significantly  reduced  as  the  sample 
size  is  increased;  the  Geometric  Approach  is  shown  to  have  a  lower  error  for  extremely 
large  datasets. 

14A  normal  matrix  is  m  x  n,  but  since  the  estimator  uses  k  samples  the  fastest  implementation  tested 
uses  m  x  n  x  k  matrices  and  dot  product  operations  to  reduce  time-consuming  memory  allocation  op¬ 
erations. 
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Processing  Time  vs  Number  of  Samples 


Figure  5.12:  TWOSTEP  (blue).  Geometric  Approach  (red),  and  TS  +  GA  (magenta): 
processing  time  vs  sample  size. 

5. 3. 4. 3  Individual  parameter  errors 

Similarly  to  the  tests  for  changing  noise  variance.  Figures  5.15  and  5.16  show  the  Ge¬ 
ometric  Approach  marginally  outperforming  TWOSTEP  in  orthogonality  error  angle 
estimation,  and  clearly  outperforming  TWOSTEP  in  sensitivity  estimation. 

5.3.5  Summary 

The  Geometric  Approach  is  an  elegant  approach  to  the  problem  and  compares  very 
well  to  TWOSTEP,  but  suffers  from  large  memory  requirements,  and  excessive  pro¬ 
cessing  time  due  to  the  calculation  of  the  Kronecker  matrix  product  for  each  sample, 
as  well  as  some  relatively  rare  convergence  issues.  The  magnetic  calibration  facility  is 
intended  to  run  at  high  sample  rates,  with  relatively  low  noise,  allowing  suitable  esti¬ 
mation  accuracy  to  be  achieved  through  the  use  of  large  datasets  processed  using  the 
TWOSTEP  algorithm. 
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Estimation  Iterations  vs  Number  of  Samples 


- 1 - TWOSTEP 

— 0 —  Geometric  Approach 
— * —  TS  +  GA  -  Only  GA  Iterations 

10° - ‘ — . . ‘ — . . ‘ — . . ‘ — . . ‘ — . 

101  io2  io3  io4  io5  io6 

N  Samples 


Figure  5.13:  Number  of  TWOSTEP  (blue).  Geometric  Approach  (red),  and  TS  +  GA 
(magenta):  iterations  required  for  convergence  vs  sample  size. 
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Figure  5.14:  TWOSTEP  (blue).  Geometric  Approach  (red),  and  TS  +  GA  (magenta): 
RMS  of  (|B  +  noise \  —  | BesfJ-mafe | )  vs  sample  size. 


Figure  5.15:  TWOSTEP  and  Geometric  Approach:  angle  errors  vs  sample  size. 
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Figure  5.16:  TWOSTEP  and  Geometric  Approach:  sensitivity  errors  vs  sample  size. 
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In  this  section  the  TWOSTEP  algorithm's  use  is  extended  to  encompass  the  calibration 
of  coils  and  sensors  within  coils.  Additionally  an  AC  calibration  is  presented,  which 
provides  sensitivity  and  orthogonality  without  estimating  bias.  In  the  first  section  (6.1), 
the  coil  to  be  calibrated  refers  to  the  relatively  large  coil  (e.g.  Ternan  coil),  into  which 
a  calibrated  reference  sensor  is  placed.  In  the  second  section  (6.2),  the  un-calibrated 
sensor  is  calibrated  by  placing  it  within  the  previously  calibrated  large  coil. 


6.1  Coil  Calibration  Estimation 


Coil  calibration  is  preferably  done  with  an  Overhauser  sensor,  hence  the  total  field 
must  always  remain  within  the  sensor's  operating  range.  This  can  be  performed  using 
earth's  field  as  a  "bias"  field  placing  the  sensor  into  its  operating  range,  and  applying 
only  signal  amplitudes  which  do  not  take  the  Overhauser  sensor  out  of  its  operating 
range.  Alternatively  a  calibrated  fluxgate  vector  magnetometer  can  be  used  as  the 
reference  sensor,  at  the  expense  of  a  lesser  precision  of  results. 


6.1.1  Model 

We  begin  by  defining  the  coil  model  as  follows: 


( I  +  D)SIk  +  Ae(be  +  4)  =  dTAsHk  +  4  (6.1) 


Page  101 


6.1  Coil  Calibration  Estimation 


where  Ik  is  the  kth  reading  of  the  current  driving  the  coils15,  and  S  is  a  diagonal  matrix 
representing  the  coil  factors.  D  is  a  fully  populated  symmetric  matrix  representing  coil 
factor  and  orthogonality  errors  when  added  to  the  identity  matrix  I.  Hk  is  the  refer¬ 
ence  magnetometer  reading,  be  is  the  bias  magnetic  field  associated  with  earth's  field, 
eek  is  the  earth's  geo-magnetic  noise,  and  e'k  is  reference  sensor  (e.g.  Overhauser)  noise. 
Ae  is  an  attitude  correction  matrix  bringing  earth's  field  into  the  coil's  reference  frame, 
while  As  and  Q  are  an  attitude  correction  for  the  reference  sensor  and  a  rotational  er¬ 
ror  within  the  reference  sensor.  These  terms  serve  no  part  in  the  estimation,  but  are 
included  to  show  similarity  to  the  TWOSTEP  method.  It  is  assumed  the  coil  produces 
no  magnetic  field  when  there  is  no  current  flowing  in  it,  removing  the  need  for  a  bias 
variable  associated  with  it. 

6.1.2  Manipulation  to  enable  use  of  TWOSTEP 

We  combine  S  and  hr  and  combine  Ae  with  be  and  ek  as  shown  in  (6.2),  to  obtain  Bk/  b 
and  e,  shown  in  (6.3). 


Bk  =  SIk 

b  =  -Aebe  (6.2) 

ek  =  -A‘4  +  4 

(7  +  D)Bk  -b-£k  =  dT AsHk  (6.3) 

Rearranging  (6.3),  gives  us  (6.4),  which  is  the  same  as  (5.1),  enabling  the  use  of  the 
Alonso  and  Shuster  TWOSTEP  model  re-derived  in  Section  5.1. 

Bk  =  (7  +  D)-\$TAsHk  +  b  +  £k)  (6.4) 

b,  which  estimates  earth's  field  in  the  coil's  reference  frame  can  be  discarded  after  the 
estimation  is  complete,  or  used  as  a  record  of  earth's  field  with  respect  to  the  coils. 

15Actual  measurement  is  the  voltage  across  a  resistor,  which  is  then  converted  to  current  7,  with  tem¬ 
perature  compensation. 
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6.1.3  Driving  the  coil 

When  driving  the  coil  to  stimulate  the  sensor,  we  are  only  truly  aware  of  a  desired  H 
and  estimated  b  and  D.  If  we  assume  As  and  #  to  be  known  and  to  be  identity  matrices 
then  (6.4)  can  be  used  for  stimulation  purposes. 

The  driving  current  J/c  is  found  as  7/f  =  S_1B/C,  which  can  be  used  directly  as  the  de¬ 
sired  input  reference  for  whenever  the  system  is  being  driven  by  a  PID  (Proportional 
Integration  Differentiation)  control  system,  or  in  the  case  of  higher  frequency  work, 
driving  voltage  implied  where  =  Rlj,  and  where  V/f  is  a  vector  representing  the 
voltage  driving  the  coil,  R  is  a  diagonal  matrix  defined  by  the  coil  winding  resistances, 
and  Ik  is  a  vector  representing  the  current  flowing  through  the  coil. 

6.1.4  Implementation  issues  in  TWOSTEP 

The  initial  centred  estimate  within  TWOSTEP  requires  a  "first  guess"  of  6',  which  is 
generally  set  to  zero.  Recall  that  6  =  (cE)T  (5.12b)  and  c  =  ( I  +  D)  b  (5.6b),  i.e.  6  is  in 
terms  of  b  and  the  symmetric  matrix  D.  In  the  case  of  coil  calibration,  be  is  used  to  rep¬ 
resent  the  bias  introduced  by  earth's  field.  The  initial  calculation  of  using  (5.15a),  is 
dependent  on  B  and  b,  and  since  in  the  typical  sensor  calibration  case  \B\  >>  |b|,  set¬ 
ting  b  =  0  will  have  little  impact  on  the  result.  For  the  typical  coil  excitation  amounts 
of  B  =  28.5  x  1CT6  T/A  (3.7),  when  \B\  «  |b|,  as  is  the  case  with  coil  calibrations,  the 
x,  y,  z  components  of  the  earth's  magnetic  field  at  DC  are  superimposed  on  the  rela¬ 
tively  weak  AC  field  generated  by  the  coil  (e.g.  cos (9)  Ternan  coil),  and  the  errors  are 
no  longer  negligible,  hence  b  should  be  as  close  to  the  real  value  as  possible.  In  this 
context,  an  un-calibrated  fluxgate  sensor  should  provide  a  suitable  first  guess  for  (b)e. 
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6.2  Sensor  in  Coil  Estimation 


The  calibration  of  sensors  in  a  coil  is  preferably  performed  with  a  reference  sensor 
co-located  with  the  sensor  under  test.  Earth's  field  cancelation  could  otherwise  be 
achieved  by  applying  an  offset  field  be,  estimated  in  Section  6.1. 

6.2.1  Model 

We  begin  by  defining  the  sensor  model  as  follows: 

Bt  =  (I+  D)-\t)TAs(A‘(H‘k  +  4')  +  HQ  +  b  +  4)  (6.5) 

where  B/(  is  the  kth  sensor  reading,  Hjr  is  earth's  magnetic  field  at  that  reading,  Hjf  is 
the  coil-generated  magnetic  field  at  that  reading,  b  is  sensor  bias,  ef  and  e'k  are  the 
measurement  noise  associated  with  earth's  field  and  the  coil  system  respectively.  D  is 
a  fully  populated  symmetric  matrix  representing  sensitivity  and  orthogonality  errors 
when  added  to  the  identity  matrix  I.  Ae  is  an  attitude  correction  matrix  (rotation),  to 
put  Hk  in  the  same  reference  frame  as  Hjf,  while  As  is  an  attitude  correction  matrix 
(rotation)  to  put  AeHek  and  Hj.  in  the  same  reference  frame  as  B/(.  i9  is  a  rotational  error 
within  the  sensor,  hence  is  another  rotation  matrix. 

We  zero  out  earth's  field,  either  by  using  a  reference  sensor,  or  the  estimated  be  from 
Section  6.1,  removing  Hjr. 

Bk  =  (I  +  D)-\«tAs(A‘4  +  HQ  +  b  +  4)  (6.6) 

Assuming  isotropic  noise  and  setting  ek  =  AsAeek,  and  =  aek  +  ek  gives  (6.7), 
which  can  be  solved  using  TWOSTEP. 

Bt=(I  +  D)-\t>TAsHct  +  b  +  ek)  (6.7) 

Hjr  and  Ae  are  expected  to  be  known  during  the  calibration,  but  only  the  |  Hk  \  needed 
for  TWOSTEP,  and  &  is  not  known  or  estimated.  As  is  likely  to  be  an  identity  matrix, 
depending  how  the  sensor  is  placed  within  the  coil. 
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Since  the  noise  has  been  merged  into  e^,  its  mean  and  variance  should  be  calculated, 
with  eek  ~  jV( 0,E|')  and  e'k  ~  A/"(0,E[.),  showing  e*.  ~  A/"(0,  E^),  where  E(e£)  = 
E{dTAsAee%)  =  0  and  Ejr.  =  Ef 


e(£,)  =  e(4)  +  e(4) 

=  0  +  0  (6.8) 

=  0 


E k  =  E(e2k)  -  E(ek) 

=  E(4) 

-i  N 

=  ^E(4  +  4)2 

iv  fc=i 

i  N 

=  n  E4  +42 +244 

fc=i 


(6.9) 


6.2.2  Bias  issues 

Due  to  slow  varying  geomagnetic  noise,  perfect  zeroing  of  earth's  field  cannot  be  as¬ 
sured,  hence  the  remaining  Hk  will  "bleed  through"  the  estimation  to  b,  introducing  an 
additional  bias  error.  This  bias  error  can  be  eliminated  if  a  calibrated  reference  sensor 
is  collocated  with  the  device  under  test,  providing  true  total  field  measurements. 
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7.1  Introduction 


This  section  presents  the  results  of  a  series  of  experiments  which  were  conducted  in 
December  2012  and  March  2013  to  perform  the  following  list  of  tasks: 

•  Magnetometer  spin  calibration  -  Calibration  of  Bartington  sensors  in  earth's  field. 

•  Calibration  for  various  sample  sizes  -  Evaluate  estimator,  and  select  suitable  cal¬ 
ibration  timeframe. 

•  Calibration  with  filters  -  Investigate  the  impact  of  filtering  on  calibration. 

•  Calibrated  Excitation  -  Excite  the  coil  using  a  calibration  matrix. 

•  Transfer  function  -  Measure  the  frequency  response  of  the  Ternan  coil. 

•  Uniformity  measurements  -  Investigate  the  usable  uniform  volume  in  the  Ternan 
coil. 

7.1.1  Experimental  setup 

The  overall  experimental  setup  is  shown  in  Figure  7.1. 

A  laptop  running  custom  Lab  VIEW  control  and  acquisition  software  is  interfaced  to 
the  system  via  two  separate  National  Instruments  acquisition  chasses.  The  first  is  a 
cDAQ  9188  chassis  containing  one  16-bit  analog  output  module  to  control  the  system, 
and  two  24-bit  differential  analog  input  modules  to  sample  the  voltage  signals  from  a 


Page  107 


7.1  Introduction 


National  Instruments  NI-ENET 
NI-9205  16  Bit  Analog  Input 


Control  Program 

LabVIEW  2012  Developer  Suite 

Runing  on  Laptop 


Bipolar  Power  supplies 
3x  Kepco  BOP  100-4ML 


■4T 


10  x  Bartington  Mag-03  MS 


1  x  AMS308i  National  Instruments  cDAQ  9188 


N I -9264  16  Bit  Analog  Output 
NI-9239  24  Bit  Analog  Input 
N 1-9229  24  Bit  Analog  Input 


Figure  7.1:  Experimental  setup  used  throughout  the  experiments. 


reference  sensor  and  reference  lO  resistors.  The  second  chassis  is  a  NI-ENET  single 
module  chassis,  containing  a  16-bit  32-channel  analog  input  module,  which  is  used 
to  sample  the  voltage  signals  from  an  array  of  10  sensors  during  the  uniformity  ex¬ 
periment.  All  analog  input  and  output  modules  are  sampled  at  6.25  kHz,  with  the 
acquisition  starts  syncronised  by  a  hardware  trigger. 

Three  Kepco  bipolar  BOP-100  amplifiers  were  used.  The  two  new  BOP-100  ML  were 
used  to  drive  the  y-axis  and  z-axis  (Ternan  design),  while  an  older  BOP-100  M  was 
used  to  drive  the  x-axis  (Solenoidal  design).  The  amplifiers  were  each  connected  to  the 
coils  in  series  with  a  IQ  50  W  reference  resistor  respectively  The  reference  resistors 
were  always  on  the  low-side  of  the  amplifiers  for  safety  reasons. 

The  reference  sensor  was  powered  by  a  Bartington  signal  conditioning  unit  (SCU),  but 
the  unconditioned  output  was  used  in  all  experiments  except  the  "Calibration  with 
filters"  dataset.  The  reference  sensor  was  always  located  at  the  coil  centre  (0,0,0),  the 
sensor  is  a  Bartington  MAG3  -  SM,  with  serial  number  "386".  The  array  of  sensors  was 
mounted  in  a  purpose  built  rig,  placing  the  sensors  in  a  line  parallel  to  the  y-axis  and 
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at  a  height  z  =  0  (Table  7.1)  and  the  array  was  moved  along  the  x-axis  using  a  trolley 
designed  for  the  purpose.  The  trolley  and  sensor  array  assembly  can  be  seen  in  Figure 
7.2. 

Since  the  earth's  magnetic  field  was  of  no  interest  in  these  measurements,  data  was 
pre-processed  by  removing  the  mean  from  each  dataset,  to  reduce  the  probability  of 
the  TWOSTEP  estimator  detecting  a  local  maxima  instead  of  a  global  maxima. 


Figure  7.2:  Array  of  Bartington  MAG3-SM  sensors  used  to  measure  field  uniformity. 

7.1.2  Stimulation  signals 

For  accurate  TWOSTEP  estimate  results,  the  reference  sensor  needs  to  be  stimulated 
by  a  signal  of  approximately  constant  "Total  field"  amplitude,  with  the  field  direction 
uniformly  distributed  on  the  surface  of  a  sphere.  Failure  to  sufficiently  spread  the  field 
direction  can  cause  the  TWOSTEP  estimator  to  generate  a  D  +  I  calibration  matrix  con¬ 
taining  imaginary  components,  which  is  incompatible  with  the  known  magnetic  field 
properties.  Hence  the  field-generating  coils  were  stimulated  using  a  signal  starting 
from  pointing  up  along  the  z-axis,  spiralling  to  point  downward,  and  then  continuing 
up  to  the  original  direction  as  illustrated  in  Figures  7.3  and  7.4. 

The  MATLAB  code  for  generating  a  200  data-point  test  signal  is  shown  below: 

N=200; 
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Sensor  Name 

Sensor  S/N 

X  Position 

Y  Position  (m) 

Z  Position  (m) 

SI 

388 

Measured 

-0.43 

0 

S2 

418 

Measured 

-0.33 

0 

S3 

392 

Measured 

-0.23 

0 

S4 

389 

Measured 

-0.13 

0 

S5 

419 

Measured 

-0.045 

0 

S6 

391 

Measured 

0.045 

0 

S7 

420 

Measured 

0.13 

0 

S8 

421 

Measured 

0.23 

0 

S9 

390 

Measured 

0.33 

0 

S10 

385 

Measured 

0.43 

0 

Table  7.1:  Sensor  array  serial  numbers  (S/N)  and  placement. 


theta= [0 : 1/N : 1] *2*pi ; 

Ratio=20; 
phi=theta*Ratio ; 

Z=cos (theta) ; 

X=sin(theta) . *sin(phi) ; 

Y=sin(theta) . *cos(phi) ; 

Ramp-up  and  ramp-down  of  the  amplifiers  were  achieved  by  applying  a  Tukey  win¬ 
dow  to  the  whole  test  signal. 
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Calibration  signal  used  in  experiments 


Figure  7.3:  Ideal  distribution  of  data-points  forming  the  calibration  test  signal. 

7.2  Magnetometer  Spin  Calibration 

7.2.0. 1  Aim 

To  calibrate  sensors  in  a  uniform,  low-noise  environment,  to  enable  latter  use  in  coil 
experiments. 

7.2.0. 2  Method 

Magnetometers  were  taken  to  a  low-noise  magnetic  facility  located  at  the  Defence  Es¬ 
tablishment,  Orchard  Hills  in  Sydney.  Sensors  were  then  individually  mounted  on  an 
approximately  1.5  m  long  fibre-glass  pole,  the  pole  was  pointed  vertically  toward  the 
sky,  and  rotated  such  that  the  sensor  remained  at  a  constant  height.  This  procedure  was 
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Calibration  signal  used  in  experiments 


Figure  7.4:  Ideal  stimulation  signals  applied  to  the  coils. 


repeated  for  each  axis,  with  a  different  axis  pointing  skyward  in  each  iteration.  Figure 
7.5  shows  the  measured  data  superimposed  on  a  sphere  of  radius  equal  to  the  earth's 
magnetic  field  scalar  magnitude  in  nanoTesla.  Data-sets  were  concatenated  prior  to 
processing  using  the  TWOSTEP  algorithm. 

The  TWOSTEP  matrices  were  reduced  using  QL  decomposition  and  angular  informa¬ 
tion  was  recovered  from  the  lower-triangular  matrix  form.  Results  are  presented  below 
in  Tables  7.2,  7.3,  and  7.4. 

Results  proved  to  be  generally  within  the  manufacturer-claimed  specifications.  Sensor 
386  was  selected  for  use  as  the  reference  sensor  in  the  coil  measurements. 
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Calibration  data  -  S/N:386  (500  nT) 
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Figure  7.5:  Calibration  data  obtained  by  spinning  the  reference  sensor  about  it  axes  in 
a  uniform  magnetic  field. 

7.3  Calibration  Experiments 

7.3.1  Coil  calibration 

7. 3. 1.1  Aim 

To  excite  the  coil  with  sample  sets  of  increasing  length,  to  investigate  convergence  of 
calibration  factors  to  their  true  values  as  sample  set  size  is  increased. 

7.3.1. 2  Method 

The  coil  system  was  connected  as  detailed  in  Figure  7.1,  with  the  sensor  array  placed  at 
the  end  of  the  coil.  The  coil  was  allowed  to  warm  up  for  a  few  minutes  by  running  an 
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SensorS /N 

Kx 

Ky 

Kz 

159 

0.99374 

0.98745 

0.99009 

161 

0.99171 

0.99171 

0.98968 

315 

1.00197 

0.99913 

0.99942 

316 

1.00332 

0.99979 

1.00027 

317 

1.00260 

1.00210 

0.99966 

322 

0.99314 

0.99148 

0.99409 

386 

1.00009 

1.00565 

1.00448 

387 

1.00462 

1.00379 

1.00497 

388 

1.00429 

1.00558 

1.00383 

389 

1.00478 

1.00286 

1.00315 

390 

1.00560 

1.00493 

1.00306 

391 

0.99902 

0.99842 

1.00108 

392 

1.00137 

0.99966 

1.00325 

418 

1.00070 

1.00109 

1.00124 

419 

0.99998 

0.99928 

0.99985 

420 

1.00018 

1.00030 

1.00088 

421 

1.00031 

0.99919 

1.00194 

Table  7.2:  Bartington  sensor  sensitivity  errors.  Shown  to  ~  5  significant  figures  consis¬ 
tent  with  Section  5.3.4  for  large  sample  sets. 


offset  current  before  and  between  measurement  runs,  and  all  measurements  are  made 
by  the  reference  sensor  located  atx  =  y  =  z  =  0 

Samples  sets  were  taken  for  the  following  durations:  1, 2, 3, 4,  5,  6,  7,  8, 9, 10, 15, 20,  30, 
40,  50,  60,  and  90  seconds. 


7. 3. 1.3  Results  &  Discussion 

Figures  7.6  and  7.7  show  the  convergence  of  the  estimated  coil  factors  as  the  sample- 
set-size  is  increased  from  1  s  to  90  s.  For  short  calibration  time-frames  it  is  believed  the 
system  exhibited  nonlinearities  associated  with  the  large  current  gradients  required  to 
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SensorS /N 

P  (de§) 

cp  (deg) 

A  (deg) 

159 

-0.34417 

-0.09965 

0.25455 

161 

-0.01075 

-0.03184 

-0.19319 

315 

-0.78405 

0.24710 

-0.09915 

316 

-0.79654 

0.09066 

0.39563 

317 

-1.29952 

0.39933 

0.07478 

322 

-0.28263 

-0.25253 

-0.13953 

386 

-0.28920 

0.16337 

-0.01491 

387 

-0.35171 

0.19582 

-0.10789 

388 

-0.12805 

0.29536 

0.05781 

389 

-0.17107 

0.22887 

-0.01739 

390 

-0.17035 

0.25021 

0.11852 

391 

-0.22311 

0.23156 

0.08445 

392 

-0.25378 

0.30264 

-0.01593 

418 

-0.27751 

0.25580 

0.02654 

419 

-0.20675 

0.20195 

-0.01028 

420 

-0.27895 

0.25246 

0.01144 

421 

-0.23791 

0.33610 

-0.13507 

Table  7.3:  Bartington  sensor  orthogonality  errors.  Shown  to  ~  5  significant  figures 
consistent  with  Section  5.3.4  for  large  sample  sets. 


drive  the  coils  with  the  calibration  signals  in  a  compressed  time-frame.  For  larger  data¬ 
sets  the  estimator  outputs  converged  to  the  results  in  Figure  7.6,  which  were  within  1% 
of  the  theoretically  calculated  coil  factors,  validating  the  modelling  in  Section  3.1.2. 

Figure  7.8  shows  the  convergence  of  the  orthogonality  errors  as  the  sample-set-size  is 
increased.  Similarly  to  the  Coil  Factors,  the  Orthogonality  estimates  converge  as  the 
sample-set-size  is  increased,  but  not  as  rapidly  as  the  coil  factors.  The  errors  are  within 
the  expected  fabrication  errors  of  the  Ternan  coil,  and  better  estimates  are  believed  to  be 
obtainable  if  larger  sample-sizes  are  taken  by  either  increasing  the  sampling  frequency 
up  to  the  50  kFIz  maximum,  or  increasing  the  calibration  duration  to  a  few  minutes. 
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SensorS /N 

bx  (nT) 

by  (nT) 

bz  (nT) 

159 

33.57672 

22.35900 

0.41812 

161 

12.03960 

16.44380 

-4.13686 

315 

-174.50699 

-156.40457 

-178.95622 

316 

-147.12588 

-77.25312 

-126.21145 

317 

-230.59581 

-214.92621 

-306.58141 

322 

-10.28013 

5.73829 

-2.82243 

386 

-185.77491 

-219.63206 

-197.52151 

387 

-234.02717 

-254.95418 

-229.84132 

388 

-160.93132 

-158.60122 

-208.77428 

389 

-258.03338 

-236.03616 

-199.25907 

390 

-187.44494 

-227.44446 

-201.38858 

391 

-238.83838 

-255.40458 

-221.76330 

392 

-231.59662 

-232.21947 

-223.22166 

418 

-187.21094 

-207.46959 

-206.52609 

419 

-187.64405 

-231.71029 

-221.83179 

420 

-221.85912 

-217.47157 

-211.56249 

421 

-214.83496 

-207.11894 

-220.55367 

Table  7.4:  Bartington  sensor  bias  errors.  Shown  to  ~  5  significant  figures  consistent 
with  Section  5.3.4  for  large  sample  sets. 


Out  of  curiosity  a  TWOSTEP  calibration  was  performed  between  the  theoretical  drive 
voltage  signal  at  the  amplifier  output,  and  the  current  passing  through  the  coils.  It 
was  found  the  output  represented  a  matrix  with  diagonal  elements  equal  to  the  series 
resistance  of  the  coils  and  reference  resistors,  and  the  non-diagonal  elements  showing 
the  coupling  between  the  respective  coil  axes.  Similarly  to  the  coil  factors  these  results 
converged  very  rapidly,  as  shown  in  Figure  7.9. 

From  the  above  results  it  is  evident  a  10  s  calibration  converges  sufficiently  to  enable 
comparison  of  datasets  as  other  parameters  are  varied,  but  for  a  true  coil  calibration,  a 
30  to  90  s  duration  or  greater  is  recommended. 
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Figure  7.6:  Coil  factors  converging  to  a  stable  value  as  sample-set-size  is  increased. 


7.3.2  B  field  predictions 

7.3. 2.1  Aim 

To  predict  the  generated  field  based  on  the  drive  signal  (voltage  or  current)  and  the 
calibration  matrix,  using  the  stimulation  signal  described  in  Section  7.1.2  Figures  7.3 
and  7.4. 


7. 3. 2. 2  Method 

The  drive  signals  from  the  data-sets  measured  in  the  coil  calibration  experiments  were 
re-processed  using  the  calibration  matrix  (D  +  I)  derived  from:  the  90  second  data-set, 
each  individual  data-set  and  the  calibration  matrix  defined  by  the  coil  factors  theoreti¬ 
cally  derived  in  Section  3.1.2. 
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Convergence  of  coil  factor  estimates  v.s.  sample  set  duration 


Figure  7.7:  Coil  factors  normalised  with  respect  to  the  values  estimated  by  the  90  s 
dataset. 

7. 3. 2. 3  Results  &  Discussion 

Figure  7.10  shows  the  RMS  value  of  the  magnetic  field  measured  by  the  reference  sen¬ 
sor  (located  at  the  centre  of  the  Ternan  coil)  taken  away  from  the  predicted  field  based 
on  the  theoretically  derived  coil  factors,  the  calibration  matrix  derived  for  the  individ¬ 
ual  data-sets,  and  the  calibration  matrix  derived  using  the  90  second  data-set.  The  90 
second  data-set  calibration  matrix  outperforms  the  theoretically  derived  matrix,  and 
matches  the  results  for  the  data-set  specific  matrices  at  large  duration  data-sets  and 
outperforms  them  for  shorter  time-frames.  Figures  7.11  and  7.12  show  the  error  in  the 
predicted  fields,  for  the  4  and  40  second  data-sets,  they  appear  to  indicate  the  greatest 
prediction  errors  are  where  there  are  the  greatest  magnetic  field  gradients  in  the  signal. 
This  is  likely  due  to  small  phase  offsets  between  the  signals,  being  most  evident  for  fast 
changing  signals.  Note  the  large  oscillations  in  the  theoretical  curves  (green)  in  Figures 
7.11  and  7.12  ,  are  caused  by  the  theoretical  coil  factors  deviating  from  reality. 
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Figure  7.8:  Orthogonality  errors  converging  to  stable  values  as  set-size  is  increased. 
Convergence  is  not  as  rapid  as  that  of  the  coil  factors. 


Figures  7.13,  7.14  and  7.15  show  the  output  field  can  be  more  effectively  predicted 
from  the  measured  current  flowing  through  the  coil,  than  from  the  voltage  driving  the 
signal  amplifiers.  This  is  mostly  likely  due  to  small  fluctuations  of  either  the  drive 
voltage  seen  by  the  amplifier,  or  of  the  amplifier  output.  Any  voltage  fluctuations  in 
the  amplifier  were  captured  by  the  current  measurement  leading  to  a  more  accurate 
prediction  of  the  generated  magnetic  field. 

This  means  for  high  precision  magnetic  field  generation,  a  feedback  system  will  still  be 
necessary  to  maintain  a  desired  magnitude,  and  compensate  for  minor  fluctuations. 

Sudden  jumps  in  the  magnetic  field  predictions  with  respect  to  voltage  shown  in  Figure 
7.15  have  been  tracked  down  to  originating  from  the  x-axis  coil  prediction.  Figure  7.16 
shows  a  sudden  jump  in  the  measured  current  which  coincides  with  the  jumps  in  the 
predicted  field.  This  is  believed  to  be  an  effect  of  the  amplifier  not  being  optimised  to 
drive  high  inductance  loads. 
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Estimated  resistance  of  individual  coils  v.s.  sample  set  duration 
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Figure  7.9:  Estimated  system  resistance  viewed  from  the  amplifier  output  as  sample- 
set-size  is  increased. 


The  remaining  noise  appears  to  be  principally  50  Hz  noise  introduced  into  the  system 
by  un-shielded  wires  driving  the  amplifiers.  The  noise  floor  within  the  facility  was  17.5 
nT  peak  to  peak. 


7.3.3  Analog  filter  impacts 

7.3.3. 1  Aim 

To  excite  the  coil  with  and  without  filters,  to  investigate  changes  in  calibration  accuracy 
if  noise  is  reduced. 

7. 3. 3.2  Method 

The  coil  system  was  connected  as  detailed  in  Figure  7.1,  with  the  sensor  array  placed 
at  the  end  of  the  coil.  The  coil  was  allowed  to  warm  up  for  a  few  minutes  by  running 
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Estimated  B<7enire|  error  based  on  sensed  current  v.s.  Sample  set  duration 


Figure  7.10:  RMS  error  of  field  prediction  using  theoretical  and  estimated  calibration 
matrices. 

an  offset  current  before  and  between  measurements.  Measurements  were  taken  at  the 
Bartington  SCU-1  filter  settings  shown  in  Table  7.5. 

Note  that  in  this  and  subsequent  measurements  sets,  the  sense  resistors  for  the  x-axis 
and  y-axis  were  shorted  due  to  a  loose  connection,  the  resulting  coupling  changes  the 
resultant  coil  factors  with  respect  to  the  drive  voltage,  since  the  system  resistance  is 
reduced  by  ~  0.50,  and  there  is  an  increase  in  coupling  between  the  two  axes.  Hence 
the  subsequent  figures  are  referenced  to  the  system  drive  voltage,  and  whilst  the  coil 
factors  and  couplings  are  slightly  different  to  those  measured  prior  to  the  reference 
resistors  shorting,  the  results  are  not  invalidated. 

7. 3. 3. 3  Results  &  Discussion 

Figures  7.17  and  7.18  show  the  estimated  orthogonality  errors  and  coil  factors  for  each 
of  the  filter  settings  respectively.  It  is  evident  the  bandpass  filters  cause  the  estimator 
to  generate  spurious  results,  this  is  likely  caused  by  the  oscillations  introduced  by  the 
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Figure  7.11:  \Bcentre\  —  using  theoretical  and  estimated  calibration  matrices. 

4  second  calibration. 


Description 

LPF  (Hz) 

HPF  (Hz) 

APF 

N/A 

N/A 

LPF1 

1000 

0 

LPF2 

100 

0 

LPF3 

10 

0 

LPF4 

1 

0 

BPF1 

10 

0.01 

BPF2 

100 

0.01 

BPF3 

1,000 

0.01 

BPF4 

10,000 

0.01 

Table  7.5:  Bartington  SCU-1  filter  settings  (3dB)  used  in  experiment. 
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Figure  7.12:  \Bcentre\  —  \Bestimated\  using  theoretical  and  estimated  calibration  matrices. 
40  second  calibration. 

bandpass  filters  whenever  there  is  a  large  gradient  within  the  field,  as  shown  in  Fig¬ 
ure  7.19.  The  LPF4  data-set  generated  incorrect  results  due  to  the  filter  removing  the 
majority  information  present  in  the  stimulating  field. 

Focusing  in  on  data-sets  APF,  and  LPF1-3,  Figures  7.20  and  7.21  show  the  filtering  does 
not  have  a  major  impact  on  the  estimation  of  coil  factors,  but  the  angles  are  affected  as 
the  low-pass  filter  corner  frequency  is  lowered.  LPF3,  the  only  filter  which  eliminates 
50  FIz  noise  also  affects  the  orthogonality  error  estimation,  increasing  the  estimated 
errors. 

From  Figure  7.22,  it  appears  filters  LPF1-3  all  track  the  APF  extremely  well,  but  if  ob¬ 
served  close  up  as  shown  in  Figures  7.23  and  7.24,  the  filtering  causes  significant  di¬ 
versions  from  the  true  signal  at  inflection  points,  and  generally  a  small  delay  in  the 
signal. 
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Estimated  \Bcentre\  error  v.s.  Sample  set  duration 
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Figure  7.13:  RMS  error  of  field  prediction  using  estimated  voltage  and  current  refer¬ 
enced  calibration  matrices. 
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Figure  7.14:  |Bce„tre|  —  estimated  using  voltage  and  current  referenced  calibration  ma¬ 
trices.  4  second  calibration. 


Page  125 


7.3  Calibration  Experiments 


fH 

s 


)5 

CQ 


e 

| 


m 


Estimated  \Bcentre\  error 


100 


50 


-50 


-100 


-150 


Current  Reference 
Voltage  Reference 


10  15  20 

Time  (s) 


25  30  35 


Figure  7.15:  \Bcentre\  —  ^estimated]  using  voltage  and  current  referenced  calibration  ma¬ 
trices.  40  second  calibration. 


Possible  jumps  in  amplifier  output 


Figure  7.16:  Unexpected  step  increase  in  the  current  flowing  in  the  x-axis  coil. 


Page  126 


Chapter  7 


Experimental  Results 


Estimated  orthogonality  error  for  filtered  signals 


Figure  7.17:  Estimated  orthogonality  errors  within  the  Ternan  coil  system,  using  a  va¬ 
riety  of  analog  filters. 
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Figure  7.18:  Estimated  coil  factors  within  the  Ternan  coil  system,  using  a  variety  of 


analog  filters. 


Page  128 


Chapter  7 


Experimental  Results 


Figure  7.19:  Recorded  \Bcentre\  data,  using  a  variety  of  analog  filters. 
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Estimated  orthogonality  error  for  filtered  signals 


1.5 


Figure  7.20:  Estimated  orthogonality  errors  within  the  Ternan  coil  system,  using  an 
all-pass-filter  and  some  low-pass-filters. 
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Figure  7.21:  Estimated  coil  factors  within  the  Ternan  coil  system,  using  an  all-pass-filter 
and  some  low-pass-filters. 


Estimated  coil  factor  for  filtered  signals 


Sensed  |B|  vs.  filter  settings 


Figure  7.22:  Recorded  \B centre]  data,  using  an  all-pass-filter  and  some  low-pass-filters. 
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Sensed  |B|  vs.  filter  settings 


Figure  7.23:  Filter-induced  deviation  from  the  unfiltered  (APF)  magnetic  signal. 


Sensed  |B|  vs.  filter  settings 


Figure  7.24:  Filter-induced  temporal  delay  of  signals. 
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From  the  above  data  it  is  evident  filtering  has  the  potential  to  adversely  affect  the  es¬ 
timator,  but  when  the  cut-off  frequency  is  lowered  sufficiently  to  remove  the  predom¬ 
inant  noise  source  (50  Hz  interference)  the  calibration  signals  are  also  affected.  An 
extension  of  the  calibration  signal  duration  should  lower  the  calibration  signal  com¬ 
ponent  frequencies  sufficiently  to  be  unaffected  by  the  filters,  but  for  such  signal  du¬ 
rations  the  estimators  will  have  already  converged  to  stable  values  due  to  the  sample- 
set-size. 

While  it  sounds  counter-intuitive  it  would  appear  the  estimator  performs  better  with¬ 
out  any  filtering,  or  only  sufficient  filtering  to  remove  signals  above  the  Nyquist  fre¬ 
quency  of  the  sampling  system. 
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7.4  Coil  Property  Measurements 

7.4.1  Coil  frequency  response 

7.4.1. 1  Aim 

To  measure  the  generated  magnetic  field  with  respect  to  the  drive  voltage  as  the  fre¬ 
quency  is  changed. 

7.4.1. 2  Method 

The  coil  system  was  connected  as  detailed  in  Figure  7.1,  with  the  sensor  array  placed 
at  the  end  of  the  coil.  The  coil  was  allowed  to  warm  up  for  a  few  minutes  by  run¬ 
ning  an  offset  current  before  and  between  measurements.  The  frequency  was  changed 
from  1  to  90  Hz,  stimulating  each  axis  coil  sequentially,  generating  independent  sets  of 
measurements  for  each  coil. 

7. 4. 1.3  Results  &  Discussion 

Figure  7.25  shows  the  coil  factor  with  respect  to  the  drive  voltage  decreases  as  the  drive 
frequency  is  increased.  Figure  7.26  shows  the  3  dB  point  of  the  system  is  80  Hz.  The 
attenuation  observed  in  Figures  7.25  and  7.26  is  due  to  the  amplifier's  inability  to  drive 
the  inductive  load  at  high  frequencies. 

Figure  7.27  shows  the  phase  response  of  the  system  between  the  drive  voltage  at  the 
amplifier  input,  and  the  reference  sensor  as  the  frequency  is  increased.  The  phase 
change  is  principally  attributable  to  a  ~  8  ms  delay  introduced  by  the  amplifier. 
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Coil  Factor  v.s.  frequency 


Figure  7.25:  Coil  factors  for  Ternan  coil  system  as  frequency  is  increased  from  1  FIz  to 
90  FIz.  Coil  factors  are  with  respect  to  the  amplifier  drive  signal. 

7.4.2  Uniformity 

7.4.2. 1  Aim 

To  measure  the  uniformity  of  the  magnetic  field  generated  by  the  magnetic  coil  system 
(Coil  length  of  8  m,  and  diameter  of  2  m). 

7.4.2. 2  Method 

The  coil  system  was  connected  as  detailed  in  Figure  7.1,  with  the  sensor  array  placed 
at  the  centre  of  the  coil.  The  coil  was  allowed  to  warm  up  for  a  few  minutes  by  running 
an  offset  current  before  and  between  measurements.  The  frequency  of  the  drive  signal 
was  set  to  3  FIz,  stimulating  each  axis  coil  sequentially,  generating  independent  sets  of 
measurements  for  each  coil.  After  each  measurement  the  array  of  sensors  was  moved 
a  small  distance  (~  5  cm)  along  the  x-axis  toward  the  end  of  the  coil  (shown  in  Figure 
7.28).  The  reference  sensor  remained  at  the  centre  of  the  coil,  to  enable  the  relative 
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[Bcentre]  Attenuation  v.s.  frequency 


Figure  7.26:  Attenuation  of  generated  magnetic  field  for  the  Ternan  coil  system  as  fre¬ 
quency  is  increased  from  1  Hz  to  90  Hz. 


variation  in  field  amplitude  to  be  measured  as  the  array  of  sensors  is  moved  along  the 
x-axis. 


7.4.2. 3  Results  &  Discussion 

Figure  7.29  shows  the  x-axis  coil  has  a  large  ripple  on  it  correlated  to  the  spacing  of  the 
coil  segments  forming  the  solenoid,  with  the  sensors  closest  to  the  boundary  showing 
the  ripple  most  strongly.  Figures  7.30  and  7.31  show  the  y-axis  and  z-axis  coils  exhibit 
smooth  uniformity  transitions  as  the  sensors  are  moved  along  the  x-axis,  with  a  large 
drop-off  at  ~  3  m  from  the  coil  centre. 

Figure  7.32  shows  the  magnitude  of  the  magnetic  field  generated  by  each  coil,  for  situa¬ 
tions  where  only  magnetic  field  amplitude  uniformity  is  necessary,  this  graph  indicates 
a  cylindrical  object  up  to  2.5  m  long  and  50  cm  in  diameter  can  be  placed  into  a  0.5% 
uniformity  field,  assuming  symmetric  field  uniformity  along  the  x-axis. 
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|Bcentre|  Phase  v.s.  frequency 


Figure  7.27:  Phase  delays  in  the  Ternan  coil  system  as  frequency  is  increased  from  1  Hz 
to  90  Hz.  Delays  are  with  respect  to  the  amplifier  drive  signal. 


Figure  7.28:  Trolley  and  array  of  sensors  at  end  of  coil. 
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|BTotai|/lBCentre|  for  X-axis  (Solenoid) 


Figure  7.29:  Total  Field  strength  for  3FIz  Bx  signals,  sensors  placed  parallel  to  y-axis 
(Table  7.1)  and  moved  along  x-axis. 


Figure  7.33  shows  a  much  more  stringent  definition  of  magnetic  field  uniformity,  where 
the  variations  from  the  desired  signal  components  are  added  vectorially,  as  in  equation 
(3.12).  In  this  case  it  appears  only  a  ~  1%  field  uniformity  can  be  obtained. 

Figure  7.34  shows  the  Bt,  By  and  Bz  fields  generated  when  only  stimulating  the  x-axis 
coil.  The  By  generated  appears  much  closer  to  zero  than  the  B;  generated,  indicating 
the  non-uniformity  is  likely  caused  by  a  slight  error  in  the  sensor  alignment  with  the 
coils,  hence  the  ~  1%  uniformity  claimed  above  is  likely  an  overestimate. 

Fig.  7.35  and  7.36  show  the  Br,  B „  and  B-  fields  generated  when  only  stimulating  the 
y-axis  and  z-axis  coils  respectively.  Note  that  the  By  and  Bz  1%  boundary  agree  with 
the  shape  generated  by  computer  modelling  in  Section  ??.  The  slight  asymmetry  along 
the  y-axis  in  Fig.  7.35  (bottom)  is  believed  to  be  due  to  a  minor  error  in  the  placement 
of  the  sensors  on  the  trolley. 
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|Brotaz|/|Bcentre|  for  Y-axis  (Ternan) 


Figure  7.30:  Total  Field  strength  for  3FIz  Bv  signals,  sensors  placed  parallel  to  y-axis 
(Table  7.1)  and  moved  along  x-axis. 

The  Ternan  coil  has  been  demonstrated  to  generate  a  large  long  volume  of  uniformity, 
and  will  be  suitable  for  a  variety  of  ranging,  calibration  and  measurement  experiments 
run  by  DSTO. 
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|BTotaz|/|Bcentre|  for  Z-axis  (Tertian) 


Figure  7.31:  Total  Field  strength  for  3FIz  Bz  signals,  sensors  placed  parallel  to  y-axis 
(Table  7.1)  and  moved  along  x-axis. 
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2-axis  (axial)  distance  from  coil  centre  (m) 


|Brota;|/|Bcentre|  for  Z-axis  (Ternan) 


x-axis  (axial)  distance  from  coil  centre  (m) 


Figure  7.32:  Total  Field  contour  plot  for  3  FIz  BX/  By  and  B2  coil  signals.  Sensors  were 
placed  parallel  to  y-axis  (Table  7.1)  and  moved  along  x-axis. 
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Figure  7.33:  Total  Field  error  contour  plot  for  3  FIz  BA/  B „  and  B;  coil  signals.  Sensors 
were  placed  parallel  to  y-axis  (Table  7.1)  and  moved  along  x-axis. 
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Figure  7.34:  Bx,  Bt/  and  B2  contour  plot  for  3  Hz  Bt  coil  signals.  Sensors  were  placed 
parallel  to  y-axis  (Table  7.1)  and  moved  along  x-axis. 
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Figure  7.35:  BZ/  By  and  Bz  contour  plot  for  3  FIz  By  coil  signals.  Sensors  were  placed 
parallel  to  y-axis  (Table  7.1)  and  moved  along  x-axis. 
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Figure  7.36:  BX/  Bv  and  B2  contour  plot  for  3  FIz  B2  coil  signals.  Sensors  were  placed 
parallel  to  y-axis  (Table  7.1)  and  moved  along  x-axis. 
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This  thesis  investigated  existing  magnetometer  calibration  estimators,  implemented 
them  in  MATLAB  and  bench-marked  their  performance  to  select  the  most  computa¬ 
tionally  fast  algorithm.The  magnetometer  calibration  models  were  extended  to  enable 
their  use  for  the  calibration  of  triaxial  magnetic  coil  systems. 


The  best  performing  estimator  was  field  tested  against  both  magnetometers,  and  mag¬ 
netic  coils,  in  both  cases  converging  to  a  fixed  solution  using  relatively  small  datasets. 
Estimated  calibration  matrices  were  shown  to  improve  the  prediction  of  generated 
magnetic  fields  when  compared  to  theoretically  derived  coil  factors.  Analog  filtering 
was  used  to  reduce  interfering  noise  sources,  but  was  found  to  introduce  phase  delays 
and  distortions  into  the  signals,  degrading  rather  than  enhancing  results. 


Additionally  this  thesis  analysed  and  modelled  the  design  of  the  previously  undoc¬ 
umented  Ternan  coil,  developing  an  improved  magnetic  coil  design  in  the  process, 
named  the  "ELFcage".  Magnetic  field  uniformity  measurements  were  conducted  on 
the  Ternan  coil,  agreeing  with  numerical  modelling  of  the  coil  system.  The  magnetic 
uniformity  was  found  to  be  sufficiently  large  to  house  most  equipment  being  mag¬ 
netically  tested  by  DSTO.  The  impact  of  high  permeability  shielding  was  also  briefly 
investigated,  with  modelling  showing  the  magnetic  field  uniformity  can  be  increased 
by  selecting  the  appropriate  shielding  geometry. 


8.1  Recommendations 


Given  the  DSTO  magnetic  test  facility  is  already  built  and  available  for  use  it  is  rec¬ 
ommended  that  the  existing  Ternan  cos(0)  coil  configuration  be  kept,  and  be  used  for 
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characterisation  of  ranging  system  sensors,  and  measurement  of  magnetic  signatures 
of  bulky  equipment.  It  is  also  recommended  a  small  2-axis  ELFcage  coil  be  built  for 
y-axis  and  z-axis  field  stimulation  of  small  sensors  within  a  shielded  laboratory  envi¬ 
ronment,  for  use  in  the  design  and  development  of  future  ranging  systems,  and  that 
an  additional  removable  solenoid  be  built  to  stimulate  the  x-axis  field. 
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A.l  Moments  of  a  Gaussian  random  variable 


First,  second,  third  and  fourth  order  moments  of  a  Gaussian  RV  are  listed  below,  as 
presented  in  problem  sheet  4  of  White  (2012).  Supposing  X  ~  N(}i,  cr2) : 

The  first  order  moment  is  shown  in  (A.l). 


E{Xf}  =  E{X} 


(A.  la) 


=  ¥ 


(A.lb) 


Second  order  moments  for  differing  and  matching  indices  are  shown  in  (A.2)  and  (A.3) 
respectively. 


E{XlXj}  =  E{X}2 


(A.2a) 


¥ 


(A.2b) 


e{x;xj  =  e{x2} 

=  }i2  +  cr2 


(A.3a) 

(A.3b) 


Page  149 


A.l  Moments  of  a  Gaussian  random  variable 


Third  order  moments  for  three  differing  indices  is  shown  in  (A.4),  two  matching  in¬ 
dices  are  shown  in  (A.5),  and  all  three  matching  indices  are  shown  in  (A.6). 


E{XiXjXk}=E{X}3 

(A.4a) 

=  v3 

(A. 4b) 

E  {XjXiXj}  =  E  |x2|  E  {X} 

(A.5a) 

=  ( }l 2  +  cr2)}i 

(A.5b) 

=  }i°  +  a2}i 

(A.  5c) 

E  {X,XfX,}  =  E  {X}3 

(A.6a) 

=  (}l3  +  3  <72}l) 

(A.6b) 

Third  order  moments  for  four  differing  indices  is  shown  in  (A .7),  two  matching  indices 

are  shown  in  (A. 8),  two  pairs  of  matching  indices  are  shown 

in  (A.9),  three  matching 

indices  are  shown  in  (A. 10),  and  all  four  indices  matching  are 

shown  in  (A.ll). 

E  {XjXjXkXi}  =  E  {X}4 

(A. 7a) 

=  / 

(A. 7b) 

E  {XiXiXjXt}  =  E  {X2}  E  {X}2 

(A. 8a) 

=  (}l2  +  cr2)ji2 

(A. 8b) 

=  }lA  +  d2jl2 

(A.8c) 

E  {XiXiXjXj}  =  E  {X2} 

(A.9a) 

=  {/  +  ^2)2 

(A.9b) 

=  ji4  +  2  a2}i2  +  u4 

(A.9c) 
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E  {XjXjXjXy}  =  E  {X}3  E  {X} 

(A. 10a) 

=  (}i3  +  3<T2pi)}l 

(A.  10b) 

=  ^(4  +  3(72//2 

(A.IOc) 

E  {XiXiXiXi}  =  E  {X4}  (Alla) 

=  }iA  +  6  cr2}i2  +  3(7 4  (A.  lib) 


A. 1.1  Fourth  moment  of  noise  amplitude 

Equation  (5.15a)  in  Section  B.1.4  contains  the  estimate  of  the  fourth  power  of  noise  sam¬ 
ple  £fc.  E(gjt)  =  }i  =  0,  and  E(£2)  —  £(£*■)  =  E*.  This  estimate  is  not  straightforward,  so 
will  be  derived  below: 


£fc)4} 

(A.  12a) 

fjOfc'ffc)2} 

(A.  12b) 

E{(4+4+4)2} 

(A.  12c) 

E  {4 + 42 + 4 + 2(44 + 44 + 44)} 

(A.12d) 

(E4  +  644  +  34)  +  2(  E  4 + 244  +  <4) 

i= 1  i=l 

(A.12e) 

3tr  (2?)  +  2tr  (z?) 

(A.12f) 

5fr  (e2) 

(A.12g) 

Assuming  noise  is  isotropic,  we  can  say  =  |e^2 1  =  | £/(3 
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EH}2=  (-tr(Zk))2  (A. 13a) 

=  (4+4  +  44  (A.  13b) 

=  4 + 4 + 4 + 2  (4  4 + 44 + 44)  (A-13c) 

=  3(4  +  4  +  4)  (A.  13d) 

=  3  tr  (A.13e) 


A. 2  Sum  of  centre  and  centred  pseudo-covariance 

Equation  (5.21)  from  5.1.7  relies  on  a  summation  over  the  product  of  centred  and  centre 
variable  to  equal  zero.  The  derivation  is  shown  below. 

N  ^  N  ^ 

E  -2(Pk-fa)(V-fl)  =  E  -2(vk-V-Hk  +  P-)(V~n)  (A.14a) 

k= 1  ak  k= 1 

N  -]  N  i 

=  E^(^-  -  ?)  -  E  tO7  -  V?  (A.  14b) 

fc— l  it— l 

N  -i  N  -] 

=  (v  -  n)  E  -  h)  -  (y  -  v)1  E  -2  (a.mc) 

^=1  ak  /c=i  ^ 

=  (A.14d) 

=  0  (A.14e) 
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A. 3  Expectations  of  centre  and  centred  variables 


/oo 

xg(x)dx.  For  E  the  probability  distribution  g(x)  is  a  delta 

-OO 

function,  so  E  {^/c}  =  /.//,.  The  same  applies  for  a2. 


{ (Vjt  -  Hfc)2}  =  E  { (v?  -  2 vkHk  +  f/?) } 


=  E  {v*}  -  E  {2w}  +  E 

=  °k  +  i4-  2nE  in}  +  vl 

=  4  +  2Vk  ~  2Vk 


=  °k 


(A. 15a) 
(A.  15b) 
(A.15c) 
(A.15d) 
(A.15e) 


A. 3.1  Centre  variables 


E  lo2 


N  -i  N  ■ i 

E  ETcT^I  ^  ^ 

=rn/  k=\ak 


k= 

a2 


(A.16a) 


(A. 16b) 


(A.16c) 


E  {/A}  =  E/c 


(A. 17a) 


Page  153 


A. 3  Expectations  of  centre  and  centred  variables 


E  M  =  E  [o2  E  \vk 

{  k= 1 

U=1  k 

N  i 

=  ^  E  “3E  W> 


=  (72 


k=l  ak 

N  i 


E  -3  n 


k=l  ak 


=  F 


E  {v  —  p}  =  E  {v}  —  E  {p} 


=  p  —  p 
=  0 


E  j  (v  —  ^)2  j  =  E  ji/2  j  —  2E  {vp}  +  E  j^2  j 
=  p2  +  a2  —  2 pE  {v}  +  p2 


A. 3. 2  Centered  variables 


E  {vk}  —  E{vk  —  v} 

=  E{vk}-E{v} 

=  pk~F 

=  Fk 


(A.18a) 

(A.  18b) 

(A.18c) 

(A.18d) 

(A.18e) 

(A.19a) 
(A.  19b) 
(A.19c) 

(A.20a) 

(A.20b) 

(A.20c) 


(A.21a) 

(A.21b) 

(A.21c) 

(A.21d) 
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E  {pk}  =  E  {}ik  -  fl}  (A.22a) 

=  E  {fik}  -  E  {fl}  (A.22b) 

=  H~V  (A.22c) 

=  (A.22d) 

E  R  -  Vk}  =  E  {vk}  -  E  {jik}  (A.23a) 

=  /A  -  /A  (A.23b) 

=  0  (A.23c) 


A. 4  Geometric  Approach  Auxiliary  Formulae 


This  section  contains  the  Auxiliary  variables  taken  directly  from  Vasconcelos  et  al. 

(2011). 

A. 4.1  Least  squares  estimate 

Scaling  an  misalignment  parameters  are  described  in  Vasconcelos  et  al.  (2011)  as: 


a  =  —  (-(4D  +  E2k2)V2  (A.24a) 

&  =  —  (-(4A  +  C2)oc2)l/l  (A.24b) 

c  =  —  ((4DA  -  B2)a2)1/2  (A.24c) 

2a\ 

tan(p)  =  ~^(2B  +  EC)(tfi)“1/2  (A.24d) 

tan (<p)  =  (BE  -  2CD)(ai)“1/2  (A.24e) 

tan(A)  =  E(—a.i^1)1^2  (A.24f) 
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Auxiliary  variables  are  defined  as: 

0!  =2 BH  +  BEI  -  2 CDI  -  4DG  +  ECH  -  E2G 
02  —  —  2 AEI  +  4 AH  -  BCI  -  2 BG  +  C2H  -  CEG 
03  =4  DIA  -  2  DGC  +  EGB  -  IB2  -  2  EH  A  +  CBH 
al  =  —  B2  +  DC2  +  4 DA  +  AE 2  -  BEC 
a2  =4 AE2J  -  E2G2  -  4 BECJ  +  2ECHG 
+  2BEIG  -  4 EHAI  -  4 DICG  -  C2H2 
+  ADAI2  +  2CBHI  -  4DG2  +  4DC2/ 

+  4 BHG  -  4 AH2  -  B2!2  -  4B2J  +  16 DAJ 
a3  =E4A  -  CBE3  +  E2C2D  -  2 B2E2 
+  8 DAE2  -  4DB2  +  16D2A 


The  error  matrix,  which  is  the  inverse  of  T/(  is  expressed  below: 


C  = 


a  0  0 

bsin(jo)  bcos(p)  0 

c(sm((p)  cos(A))  csin( A)  ccos(cp)  cos(A) 


A. 4. 2  Maximum  likelihood  estimate 


Vasconcelos  et  al.  (2011)  derived  the  log-likelihood  function  to  be: 


/(T,b)  =  E 

i=  1 


\T(hri  —  b\\  -1 
crmi 


If  we  let  Uj  =  hri  —  b,  gradient  V/ \x  is  split  into  two  sub-matrices  V/ \x 
Tk 


where  xk  = 


hi 


V/|r  =  Yj  '~rui  ®  Tui 


i= 1 


v/|„  =  t  ^At'Tu, 


■  1  (7  • 

z=l  mz 


(A.25a) 
(A.  25b) 
(A.25c) 
(A.25d) 

(A.25e) 

(A.25f) 

(A.26) 

(A.27) 

[V/|tV/|,] 

(A.28a) 
(A.  28b) 
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where  Cj  =  1  —  ||  Tui\\  1  and  0  is  the  Kronecker  product16. 
The  Hessian  V2/  \x  can  be  expressed  as: 

HT/t  HTb 
H Tb  ^ b,b 

with  sub-matrices: 


(A.29) 


Hj  j 


HT,b 

Hb,b 


i=l  u  mi 
n  —2 

i=l  u  mi 

n  2 

E 

i=l 


TiijuJ  Tt)  r/T 

—  +  CT[(UillJ)  0  Z3] 


||Tu/||3 

( Mj  0  Tui)uJttT 

||  Til;  ||  3 

TTTuiUjTTT  T 

I  Ttf,  II3  +CtTT 


+  Ct[«z  0  T  +  J3  0  TMj] 


(A.30a) 
(A. 30b) 
(A.30c) 


16A  fast  implementation  for  the  Kronecker  product  can  be  found  in  C.2.5 
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B.l  AC  calibration 


Section  6.1  alludes  to  the  coil  inherently  having  zero  bias  error  associated  with  it.  The 
only  manufacture  errors  present  are  related  to  coil  factors  and  orthogonality,  hence  b 
does  not  need  to  be  estimated  for  coil  calibration  purposes,  and  the  bias  due  to  earth's 
field  and  the  reference  sensor  would  need  to  be  removed  from  Bk  and  Hk  measure¬ 
ments.  This  could  be  achieved  by  the  use  of  a  filtering,  or  the  use  of  FFTs  as  mentioned 
in  Section  2.3.  The  following  estimator  is  effectively  a  special  case  of  TWOSTEP  in 
which  b  =  0. 


B.1.1  Model 


The  simplified  coil  calibration  model  is: 


(7  +  D)SIk  +  Aeeek  =  t 9TAsHk  +  4 


(B.l) 


Following  the  same  steps  as  shown  in  section  Section  6.1.2  gives: 


B*=(I  +  D)-1(#TAsHt  +  ek) 


(B.2) 
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Although  Hk  is  expected  to  be  known  during  the  calibration,  Ak  and  &  are  not  known 
or  estimated,  so  (B.2)  is  manipulated  to  eliminate  the  rotations  (which  are  orthogo¬ 
nal),  hence  only  magnitude  data  is  used.  Section  5.1.2  included  the  same  step,  but  for 
different  reasons. 


I  Hi:  1 2  =  ((/  +  D)Bk  -  ek)T((I  +  DjB,  -  ek)  (B.3a) 

=  Bl(l  +  D)T(I  +  D)Bt-2Bj(I  +  D)Tek  +  \et\1  (B.3b) 

If  a  scalar  sensor  were  used  to  measure  |  Hk  | ,  filtering  could  not  be  applied  to  individual 
components  of  Hk,  meaning  the  vector  sum  would  distort  the  AC  signal,  introducing 
harmonics. 

B.l. 2  Parameter  re-definition 


I  Hk  | 2  is  re-arranged  and  a  better  set  of  basis  functions  is  selected  for  improved  model 
stability  (Barker  et  al.  2007),  hence: 


\Hk\2  =  BTk  (I  +  2D  +  D2)Bk  -  2B\ (I  +  D)T ek  +  |£fc|2 

(B.4a) 

=  \Bk\2  +  BTkEBk  -  (2[(I  +  D)Bk]  ■  ek  -  \ek\2) 

(B.4b) 

=  \Bk\2  +  BkEBk  -  vk 

(B.4c) 

vk  is  defined  later  in  (B.6b).  For  ease  of  estimation  we  removed  the  non- 

dence  on  parameter  D,  defining  matrix  E  as  shown  in  (B.5). 

-linear  depen- 

E  =  2D +  D2 

(B.5) 

B.l. 3  Sensor  error  measurement  &  sensor  error  measurement  noise 

As  shown  in  (B.4),  vk  holds  all  the  ek  noise  elements  present  in  \Hk\2,  and  zk  could  be 
considered  a  scalar  measurement  of  the  sensors  errors  we  are  trying  to  estimate.  zk  and 
vk  are  defined  in  (B.6). 
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zk  =  \Bk\2  -  \Hk\2  (B.6a) 

vk  =  [(I  +  D)Bk ]  ■  ek  -  \ek\2  (B.6b) 

Using  (B.4)  the  scalar  measurement  zk  is  expressed  as  follows: 

zk  =  \Bk\2  -  \Hk\2  (B.7a) 

=  —BkEBk  +  vk  (B.7b) 

For  ease  of  derivation  and  coding  E  is  vectorised,  rearranging  to  remove  the  symmet¬ 
ric  3  x  3  matrix  E  and  to  use  6x1  vector  E  instead,  as  was  done  in  (5.9),  obtaining 
B[EBk  =  Kk E.  Where  Kk  and  E  are  defined  as  follows: 

E  —  (  Eli  E22  E33  E\2  E13  E23  )  (B.8a) 

=  (  BU  B2k  Blk  2BUB2 1-  2BlfcB3fc  2B2fcB3fc  )  <B-8b) 

Applying  the  simplification,  zk  is  expressed  in  terms  of  variables  Kk  and  E. 

Zk  =  ~BkEBk  +  vk 
=  —Kk  E  +  vk 

A  re-arrangement  of  (B.9)  represents  vk  in  terms  of  Kk  and  E  as  shown  in  (B.10). 

vk  =  zk  +  KkE  (B.10) 


(B-9a) 

(B.9b) 


B.1.4  Noise  distribution 

As  part  of  the  MLE  process  the  noise  in  the  system  needs  to  be  characterised.  As 
earlier  ek  is  assumed  to  be  white  and  Gaussian  distributed  with  zero  mean.  Since 
ek  ~  A/"(0,£fc)  then  vk  ~  N (}ik,cr2)  as  shown  in  (B.lla)  and  (B.12a). 
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Calculation  of  }ik  is  shown  below: 


H  =  E  M  (B.lla) 

=  E  {2[( I  +  D)Bk]  ■  et)  -  E  { |£i|2}  (B.llb) 

=  0-e{(^£^)2}  (B.llc) 

=  -fr(Lt)  (B.lld) 


Once  again,  calculation  of  a2  requires  the  use  of  Gaussian  random  variable  moments 
presented  in  Appendix  A.l. 


°k  =  E  { vl }  -  Fk 

=  4 E  {([(7  +  D)Bk]Tek)2}  +  2 E  {  [(/  +  D)Bk]Tek)\ek |2}  +  E  {|£fc|4} 

=  4E  { [(I  +  D)Bk]TekeTk[(I  +  D)Bk\ }  +  E  { |f,|4}  -  tr  (E,)2 
=  4 [(/  +  D)B/f]TE/f [(I  +  D)Bk\  +  5tr  (eJ)  -  3 tr  (e£) 

=  4[(I  +  D)Bk\TZk[(I  +  D)Bk]  +  2 tr  (e2) 


(B.12a) 
tr(  Efc)2 
(B.12b) 
(B.12c) 
(B.12d) 
(B.12e) 


B.l. 5  Likelihood  function 

Having  characterised  the  noise,  the  likelihood  (B.13a)  and  negative  log-likelihood  (B.14a) 
functions  are  obtained: 


N  ■ |  2 

/n|0'=e'(E')  =  n  /2”i  (B.13a) 

k- 1  v27T<7i 

N  -1 

=  ff  1  ,-(z,+K(E-ft)V2^  (6.13b) 
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/(E')  =  -ln/N|0,=E,(E') 


N  r 

E 

k=  1 

1  N 

E 

z  k=  1 


2  /o_2  ,  1  i„ 


(Zfc  +  Kl<E  ~  Jik)  /2oi  +  2  ln  (2mrfc 


-2  (Zjfc  +  KkE  -  }i) c)2  +  In  (27rcrJ 

or 


B. 1.5.1  Estimator 


The  estimator  for  the  measurement  is  found  by  setting  ^rj  =  0. 

a 


9  1^1 

9E7/  =  9F2  £  (Zfc  +  KkE  ~  HY 


+  terms  independent  of  £) 
i  N  2T<t 

=  2  E  —£  (zk  +  KkE  -  n) 
1  cri 


k= 1  uk 


N 


N  KTKi- 


ET(Zt-w)Kf+E^E  =  0 


fc=l  °"k 

N  zsT  ts  N  -i 

•■•E%?£=E  4(zk-f/fc)Xi 


fc=l 


Jt=l 


Jt=l  ak 


(B.14a) 

(B.14b) 

(B.14c) 


(B.15a) 

(B.15b) 

(B.15c) 

(B.15d) 

(B.15e) 


*  KjKk 

We  multiply  both  sides  by  the  inverse  of  ^  „ —  to  obtain  the  estimate  £*,  where  * 

fc=l 

indicates  the  value  is  an  estimate.  3  (£)  is  the  Fisher  Information  Matrix  and  is  derived 
in  Section  B.1.6 


.  £* 


-1 


E  -ji^k-  n)  KI 


k= 1  ak 


P  WP'EtT  (zt-f‘k)Kl 


(B.16a) 

(B.16b) 


B.1.6  Fisher  information  matrix 

The  Fisher  Information  Matrix  3  is  the  expectation  of  the  Hessian  matrix  -  the  square 
matrix  of  the  second  order  partial  derivatives  of  the  negative  log  likelihood  function 
/  (£).  The  fisher  matrix  is  defined  by  Poor  (1994)  as  presented  in  (5.26): 
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Using  (5.27)  the  Fisher  Information  Matrix  for  /  (£)  is  derived17,  as  shown  in  (B.17). 


3(£)=£{(^FT)} 

(B.17a) 

d2  N  1 

=  E{(^f2E  ?(zk  +  KkE  Fk)T(zk  +  KkE  fik))} 
aL  fc=i  ak 

(B.17b) 

d  N  1 

=  £{si 

3E  t=i  b 

(B.17c) 

N  -i 

=  £1E  -2KTkKk} 

(B.17d) 

k=l  ak 

N  1 

=  L  if*. 

(B.17e) 

fc= l  uk 


B.l. 7  Converting  to  D 

Estimate  E*  is  not  in  the  desired  form  D*,  but  E*  can  be  changed  back  into  matrix  form 
E*  by  using  its  definition  in  (B.8a).  To  compute  D*  we  write: 

E*  =  USUT  (B.18) 

where  U  is  orthogonal  and  S  is  a  3  x  3  diagonal  matrix  with  elements  s;,  and  both 
matrices  are  obtained  by  the  Eigen-decomposition  of  symmetric  matrix  E*. 

We  define  diagonal  matrix  W,  satisfying  S  =  2W  +  W2,  hence  using  the  quadratic 
formula  the  elements  are18: 

IV  j  =  —  1  +  y/l  +  Sj  (B.19) 

The  maximum  likelihood  estimate  for  D  can  hence  be  calculated  by: 

D*  =  UWUT  (B.20) 

17Due  to  the  differentiation  3  ( E )  is  no  longer  dependent  on  E. 

18iVj  —  —  1  —  yjl  —  Sj,  is  not  used  as  it  would  give  negative  sensitivity  values,  implying  the  sensor  was 
facing  the  wrong  way.  This  is  not  possible  since  it  would  have  been  accounted  for  by  attitude  rotations, 
and  manufacture  quality  control  is  expected  to  capture  such  major  sensitivity  errors 
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B. 1.7.1  Covariance  Matrix 

To  obtain  the  covariance  matrix  need  to  transform  3  (£)  to  3  (D),  where  D  is  defined 
with  respect  to  D  the  same  way  as  E  was  defined  in  (B.8a).  Explicit  derivation  of  this 
process  is  not  shown. 


3(D)-1 


(B.21a) 


Element  by  element  differentiation  of 


can  be  represented  as 


dEj 

dD, 


=  [MEuiD)}-1 


(B.22a) 

(B.22b) 


where 


M£D(D) 


+  2  Di 

0 

0 

2D4 

2D5 

0  N 

0 

2  +  2D2 

0 

2D  4 

0 

2  D6 

0 

0 

2  +  2  D3 

0 

2D5 

2D6 

d4 

d4 

0 

2  +  D4  +  D2 

d6 

d5 

d5 

0 

D5 

d6 

2  +  Di  +  D3 

d4 

0 

D6 

d6 

d5 

d4 

2  +  D2  +  D3  j 

(B.23) 


B.1.8  Implementation 

Begin  by  assuming  Hy,  B/c,  and  E*-  are  known. 

•  Remove  unwanted  signals  from  H/v  and  calculate  filtered  |  H/f 

•  Calculate  =  —tr( E^) 

•  Calculate  using  (B.6a). 
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•  Calculate  using  (B.8b) 

•  Set  E  =  0,  and  calculate  cr£  using  (B.12a). 

•  Calculate  Fisher  Information  Matrix  3(E*)  using  (B.17). 

•  Estimate  E*  using  (B.16). 

•  Recover  desired  values  D  as  shown  in  Sections  B.1.7  and  B.  1.7.1. 

An  implementation  of  the  code  can  be  found  in  Appendix  C.3. 

Since  the  TWOSTEP  estimator  performed  extremely  well,  it  was  decided  the  AC  Cali¬ 
bration  would  be  unnecessary,  due  to  the  large  amount  of  time  required  to  generate  a 
suitable  data-set.  Instead  the  3  EIz  frequency  bin  values  of  an  FFT  of  the  reference  sen¬ 
sor  data-sets  from  the  uniformity  measurements  in  Section  7.4.2  were  used  to  generate 
a  60  point  data-set  for  the  AC  Calibration,  the  results  of  which  are  shown  in  Figures 
B.l  and  B.2.  As  we  can  see  the  coil  factors  were  very  close  to  those  estimated  using 
TWOSTEP;  the  orthogonality  errors  are  close  to  those  estimated  by  TWOSTEP,  but  in 
one  case  of  the  wrong  sign.  The  AC  calibration  appears  to  work  as  expected,  and  may 
in  future  be  used  for  measurements  in  high  noise  environments,  where  noise  is  too 
high  for  TWOSTEP  to  handle  with  the  current  limitations  on  computer  memory. 


B.2  Amplifier  linearity 


During  the  measurement  of  the  coil  transfer  functions,  concerns  were  raised  about  the 
linearity  of  the  coil  system.  Figures  B.3,  B.4  and  B.5  were  generated  using  the  refer¬ 
ence  sensor,  and  Figures  B.6,  B.7  and  B.8  were  generated  using  the  reference  resistor 
voltages19.  It  is  apparent  that  for  each  drive  frequency  several  higher  order  harmonics 
were  generated,  with  power  ~  —  80  dB  that  of  the  fundamental. 

19Since  only  one  coil  was  being  stimulated  at  a  time,  the  frequency  distribution  measured  in  the  re¬ 
sistor  will  match  that  present  in  the  coil,  even  though  there  was  a  short  between  the  x-axis  and  y-axis 
coils. 
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Figure  B.l:  Coil  factors  estimated  by  the  AC  Calibration  algorithm  compared  to  the 
coil  factors  estimated  by  the  TWOSTEP  algorithm. 

To  identify  the  source  of  the  harmonics  the  amplifiers  were  tested  within  the  lab.  Figure 
B.9  shows  the  power  spectral  density  of  the  amplifiers  with  a  60  FIz  signal  applied. 
It  is  apparent  the  amplifiers  are  the  source  of  the  harmonics,  but  their  amplitude  is 
within  the  amplifier  specifications,  hence  whilst  there  is  a  definite  non-linearity  to  the 
amplifiers  it  is  unlikely  to  have  a  significant  impact  on  the  system. 
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Figure  B.2:  Orthogonality  errors  estimated  by  the  AC  Calibration  algorithm  compared 
to  the  orthogonality  errors  estimated  by  the  TWOSTEP  algorithm. 
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Welch  Power  Spectral  Density  Estimate  of  Bx  for  a>axis  coil 


Figure  B.3:  Measured  Welch  power  spectral  density  of  the  reference  sensor  x-axis  when 
the  coil  x-axis  is  energised. 
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Welch  Power  Spectral  Density  Estimate  of  By  for  j/-axis  coil 


Figure  B.4:  Measured  Welch  power  spectral  density  of  the  reference  sensor  y-axis  when 
the  coil  y-axis  is  energised. 
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Welch  Power  Spectral  Density  Estimate  of  B,  for  z-axis  coil 


Figure  B.5:  Measured  Welch  power  spectral  density  of  the  reference  sensor  z-axis  when 
the  coil  z-axis  is  energised. 
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Welch  Power  Spectral  Density  Estimate  of  Ix  for  s-axis  coil 


Figure  B.6:  Measured  Welch  power  spectral  density  of  current  flowing  in  the  x-axis 
coil. 
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Welch  Power  Spectral  Density  Estimate  of  Iy  for  y-axis  coil 


Figure  B.7:  Measured  Welch  power  spectral  density  of  current  flowing  in  the  y-axis 
coil. 
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Welch  Power  Spectral  Density  Estimate  of  Iz  for  2-axis  coil 


Figure  B.8:  Measured  Welch  power  spectral  density  of  current  flowing  in  the  z-axis 
coil. 
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Power  Spectral  Density  Estimate 


Figure  B.9:  Measured  power  spectral  density  of  the  Kepco  amplifiers  used  to  drive  the 
Ternan  coil  system. 
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C.l 

TWOSTEP  estimator  1 

C.1.1 

‘TWOSTEP  estimate. m’ 

function  [D_est ,  b-est ,  n,  Cov_est  ]  =TWOSTEP_estimate  (B,  H,  Sigma-noise) 
%TWOSTEP_ESTIMATE  estimation  of  magnetometer  errors. 


o, 

o 

[D-est ,  b-est ,  n,  Cov_est  ]  =TWOSTEP _e st imat e  (B,  H,  Sigma-noise ) 

o, 

o 

B  —  3xN  matrix ,  where  N  is  the  number  of  signal  samples 

Q, 

O 

H  —  Earth ' s  magnetic  field  in  Earth  frame  coordinates 

g, 

o 

Sigma-noise  —  3xN  noise  covariance  array,  where  each  column  is  a 

o, 

o 

diagonal  of  the  noise  covariance  matrix. 

g, 

o 

g, 

o 

Implemented  based  on  the  work  presented  in: 

g, 

o 

"Complete  Linear  Attitude— Independent  Magnetometer  Calibration" 

g, 

o 

by  R. Alonso  and  M.  Shuster 

g, 

o 

in  "The  Journal  of  the  Astronautical  Sciences"  Vol .  50  No. 4 

g, 

o 

October— December  2002. 

g, 

o 

g, 

o 

Sections  of  code  and  optimisation  adapted  from  the  TWOSTEP  code  written 

g, 

o 

by  J. F .Vasconcelos  for  the  validation  of  geometric  calibration  presented 

g, 

o 

in  "A  geometric  approach  to  strapdown  magnetometer  calibratio  in  sensor 

g, 

o 

frame"  by  J.F.  Vasconcelos  et  al .  in  Aerospace  and  Electronic  Systems, 

g, 

o 

IEEE  Transactions  47(2)  2011.  Many  thanks  go  to  Dr.  Jos  Vasconcelos  for 

g, 

o 

his  assistance . 

g, 

o 

g, 

o 

Slight  variation  in  the  calculation  of  the  centered  J  function,  was  used 

g, 

o 

for  clarity  in  derivations. 
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o 

o 

%  Justin  Dinale ,  DSTO  Department  of  Defence,  Australia 
%  $Revision:  1.0.0  $  $Date:  2012/11/05$ 

0.  0.0,  0.0.00.000,0000.  O  0.0.0  0  0.0,0.0.  0.0.  0.0.0.  0.0,0  0.0.0.0.0  O.  0.000,0,0,00.00.  0,0000,0,0,0.0  0,0,  0.0000,0.0.0.0.  0,0,  0.0.0. 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%  Define  stop  condition  and  perform  housekeeping 
stop-tol=le—24;  %  Stop  Condition  from  Alonso  paper 
max_n=200 ;  %  Max  number  of  iterations 
k_samples=size  (B, 2 ) ;  %  Number  of  signal  samples 


%%  TWOSTEP  Centered  estimate 
%  Set  initial  guess  for  b  and  D. 
bO=[ 0  0  0]  ' ; 

DO=zeros ( 3) ; 

%  Calculate  L 
L-k  =  [ 2  *B;  .  .  . 

-  (B.  ~2);  .  .  . 

—2  *B  (1,  :)  .  *B  (2, 

—2*B  (1,  :)  .  *B(3,  :);... 

—2*B  (2,  :)  .  *B  (3,  :)]; 

%  Calculate  scalar  measurement  z 
z.k=total-field  (B)  .  "2— total-field  (H)  .  ~ 2 ; 

%  Calculate  mean  and  variance  of  noise 
mu-k^sum  (Sigma-noise ,  1)  ;  %  Mean 

sigma-sq-k=4  *sum  (  (  (eye  (3)  +D0)  *B—bO  *ones  (1,  k.samples)  ).*... 

(Sigma-noise .  * ( (eye  (3) +D0)  *B—bO *ones  (1, k_samples ) )), 1) .. . 
+5*  (sum  ( Sigma-noise  .  ~2,  1)  )—  (sum  (Sigma-noise ,  1 )  .  ~2 )  ; 

%  Calculate  centered  sigma  squared 
s igma-Sq-bar=l / sum  (1 .  / sigma-sq-k )  ; 

%  Center  the  data 

[mu-bar,  mu -k -tilde ]  =center_data  (mu-k ,  sigma-sq-k ,  sigma-sq-bar)  ; 


Page  178 


Appendix  C 


Code 


[z_bar ,  z_k_tilde ]  =center_data  (z_k,  sigmasq.k ,  sigmasq^bar)  ; 

[L_bar,  L_k_t  i  lde  ]  =center_data  (L_k,  sigma-sq.k  ,  sigma-sq.bar)  ; 

%  Calculate  fisher  information  matrix 
[  I_fisher_tilde  ,  I_fishinv_tilde  ]  =TS_f isher_centered  ( s  igma_sq-k  ,  L_k_t  ilde)  ; 
%  Calculate  centered  estimate 
theta-0-tilde=.  .  . 

I_fishinv_tilde  *  ( L_k_ti  lde*  (  ( z_k_tilde-mu-k_tilde)  '  .  /  sigma^sq.k  ' )  )  ; 

%%  TWOSTEP  Center  correction 

theta-n=theta-0-tilde ;  %  Initiare  theta  for  first  iteration 
n=0;  %  Initialise  iteration  counter 
TS-err=Inf;  %Initial  condition  for  error. 

%  ABC  is  used  to  remove  intensive  calculations  out  of  for  loop. 

ABC=—  (  (  ( z  _k_t  i  lde— mu_k_t  ilde)  .  /  s  igma_sq-k)  *L_k_tilde  ' )  '  ; 

while  (TS-err>stop-tol  &&  n<max_n ) 

if  (n^O)  %  If  we  are  not  in  the  first  iteration 
theta-n=theta-npl ; 

end; 

%  Extract  c  and  E  components 
[ c,E]=theta-tO-C-E  (theta^n)  ; 

%  Calculate  second  derivative  of  b~2  wrt  theta 
tmp= ( (eye (3) +E) \c)  * ( (eye (3) +E)\c)  ' ; 
dbsqdtheta-p=  [2*  (  (eye  (3)  +E)  \c)  ;  —diag(tmp)  ;  .  .  . 

—2*tmp(l,2);  —2*tmp(l,3);  —  2*tmp  (2,  3)  ]  ' ; 

%Calculate  gradient  of  J 

dJdThetap-tilde=ABC  +  I-fisher-tilde*theta-n; 

dJdThetap_bar=  (—  (1/sigmasq.bar)  *  (L_bar  '  —  dbsqdtheta.p)  * .  .  . 

( z_bar  —L_bar  '  *theta_n  +c  '  *  (  (eye  (3) +E)\c)  —  mu-bar))'; 
dJdTheta=dJdThetap-tilde+dJdThetapfbar  ; 

%  Calculate  Fisher  matrix 
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[ I_fisher_bar ]  =TS_fisher_center  (sigmasq-bar ,  L_bar,  dbsqdtheta.p)  ; 

%  Update  theta 

theta.npl  =theta_n—  (I_fi  sher_tilde+I-fi  sher_bar )  \  dJdTheta; 

%  Calculate  Error 

TS-err=  (theta-npl—theta-n)  .  . 

(I_fisher_tilde  +  I-fisher_bar)  *.  .  . 

( theta.npl—theta.n)  ; 
n=n+l;  %Increase  counter 

end; 

[b-est,  D_est  ]  =theta-tO-b-D  (theta-npl )  ; 

%  Extract  covariance  matrix 
b=b.est ; 

D=D_est ; 

M.cD=[b(l)  0  0  b(2)  b  (3)  0;  ... 

0  b  (2)  0  b(l)  0  b  (3)  ;  ... 

0  0  b  (3)  0  b(l)  b  (2)  ]  ; 

M_ED=  [  2*D  (1) ,  0,  0,  2*D  (4) ,  2*D(5),  0;  ... 

0,  2*D  (2) ,  0,  2*D(4),  0,  2*D(6);  ... 

0 ,  0,  2  *D  (3)  ,  0,  2  *D  (5)  ,  2*D(6);  ... 

D  (4)  ,  D  (4 )  ,  0,  D  (1)  +D  (2)  ,  D(6),  D(5);  ... 

D  (5) ,  0,  D  (5)  ,  D  (6)  ,  D  (1)  +D  (3)  ,  D(4);  ... 

0,  D  (6)  ,  D  (6)  ,  D  (5)  ,  D(4),  D(6)+D(5)];  ... 

dbD_dcE=  (  [ (eye (3) +D) ,  M.cD;  zeros (6, 3),  M.ED ] )\eye (9) ; 

Cov_est=dbD_dcE  *  ( I_f  isher_t  ilde+I  _fisher  _bar  )  \dbD_dcE  ' ; 

%END  TWOSTEP 

function  [B_total ] =total_f ield  (B_in) 

%  Calculates  total  field  of  B_in  (3xn  matrix) . 

B_total=sum (B_in . "2,1) . "0.5; 

function  [X_bar ,  X_tilde ]  =center_data  (X,  sigma.sq-k ,  sigma.sq-bar) 
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%  Calculates  the  centered  and  center  components  of  X 
%  Center  component 

X-bar  =  sigma-sq-bar*  (X*  (1 .  / sigma-sq-k)  '); 

%  Centered  component 

X_tilde  =  X  —  X_bar *ones  (1 ,  size  (X,  2)  )  ; 

function  [  I_fisher-tilde  ,  I-fishinv-tilde  ]  =... 

TS_fisher_centered  (sigmasq,  L_tilde) 

%  Calculates  the  fisher  information  matrix  for  the  centered  estimate,  when 
%  given  variance  sigmasq  and  centered  vectors  of  L-tilde 
%  Calculate  fisher  information  matrix 
I_fisher_tilde= .  .  . 

(  (L-tilde  .  *  (ones  (size  (L_tilde,  1)  ,  1)  *  (1 .  / sigma^sq)  )  )  *L_tilde  ' )  ; 

%  Calculate  inverse 

I_fishinv_tilde=  ( I-fisher_tilde )  \eye  (9)  ; 

function  [ I_fisher_bar]  =. . . 

TS-f isher-center  ( s igma-sq-bar ,  L-bar,  dbsqdtheta-p) 

%  Calculates  center  informaiton  matris.  Used  for  readability 
I_fisher_bar=  (  (L_bar  '—dbsqdtheta-p)  '  *  (L-bar  '—dbsqdtheta-p)  )  / sigma^sq-bar ; 

function  [b,D]=  theta-tO-b-D (theta) 

%  Converts  a  value  of  theta  to  usable  physical  values 
[c,  E]  =  theta-tO-C-E  (theta)  ; 

[U,  S]  =eig  (E)  ; 

W^eye  (3)  +  (eye  (3)  +S)  .  "(0.5); 

D=U*W*U' ; 

%  Calculate  b  using  the  inversr  of  (I+D) 
b= (eye (3) +D) \c; 

function  [c,  E]  =theta-tO-C-E  (theta) 

%  Extracts  c  and  E  elements  from  theta  as  per  Alonso  paper. 
c=theta  (1:3) ; 

E= [theta  (4),  theta  (7),  theta  (8);... 
theta  (7),  theta  (5),  theta  (9);... 
theta  (8),  theta  (9),  theta  (6)]; 
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C.2  Geometric  Approach  estimator 


C.2.1  ‘ellipsoicLestimate.m’ 


function  [T,  b,  n]  =e  11  ipso  id-estimate  (h_r  ,  sigma^sq,  B_e  ) 

%ELLIPSOID_ESTIMATE  estimation  of  magnetometer  errors. 

%  [T,b,  n ]  =ellipsoid-estimate  (h_r ,  sigmasq,  B_e) 

%  h_r  —  3xN  matrix,  where  N  is  the  number  of  signal  samples 
%  sigmasq  —  IxN  noise  covariance  array,  where  each  column  is  a 
%  diagonal  of  the  noise  covariance  matrix.  Note  single  row  since  noise 
%  assumed  to  be  isotropic . 

%  B_e  total  earth  field.  IxN  vector. 

g, 

o 

g, 

o 

%  Implemented  based  on  the  work  presented  in: 

%  "Geometric  approach  to  strapdown  magnetometer  calibration 
%  in  sensor  frame " 

%  by  J.F.  Vasconcelos  et  al .  in  Aerospace  and  Electronic  Systems, 

%  IEEE  Transactions  47(2)  2011. 

o, 

o 

%  Justin  Dinale,  DSTO  Department  of  Defence,  Australia 
%  $Revision:  1.0.0  $  $Date:  2012/11/05$ 

0.0.0,  O.O.O.O.O.O.O.O.O.O.O.  O.O.O.O.O.O.O.O.O.O.O.  0.0.0.0.0,0.0,0.0.0.0.  0_0.0.0.0,0.0,0.0.0.0.  O.O.O.O.O.O.O.O.O.  O.  O.  O.  O.O.O.O.O.O.O.O.  O.  O.  0.0.0. 

oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%  Define  stop  condition  and  perform  housekeeping 
stop-tol=le—5;  %  Stop  Condition  from  Alonso  paper 
max_n=200 ;  %  Max  number  of  iterations 

k_samples=size  (h_r , 2) ;  %  Number  of  samples 

%%  Make  first— guess  estitimation  of  T  and  b 
[TO,  bO]  =ellipsoid-f irst-estimate  (h_r); 

%%  Iterative  corrections  in  estimate 

T=TO; 

b=bO; 

GC_error=Inf; 
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n=l } 

while  (abs  (GC-error)>stop-tol  &&  n<max_n) 

X=[T(:,1)  ;  T  (:,  2)  ;  T(:,3);  b]  ; 

%  If  sigma~2  is  constant,  convert  it  to  an  array  of  samples 
if  (length  (size  (sigmasq)  )  ==1) 

sigma-sq  (1,1,  : )  =sigmasq*ones  (1,  k_samples )  ; 

else 

tmp_sigma  (1,1,  :  )  =squeeze  (sigmasq)  ; 
sigma-sq=tmp_sigma ; 
clear  tmp_sigma 

end 

%%  Calculate  u_{i}  The  measurement  without  bias 
%  Calculate  temporary  vector 
u_t=h_r—  [b  (1)  *ones  (1,  k_samples )  ;  .  .  . 
b  (2)  *ones  (1,  k_samples )  ;  .  .  . 
b  (3)  *ones  (1,  k_samples )  ]  ; 

%  Put  vector  in  3D  matrix  accountign  for  samples ,  also  store  transpose 

u  (:,  1,  : )  =u_t  ; 
u -P  (1,  : ,  :)=u-t; 

%%  Calculate  regularly  used  elements 

%  Calcualte  Tu 

Tu  ( : , 1,  : )  =T * squeeze (u) ; 

%  Calculate  magnitude  of  Tu  —>  |Tu| 

ITuI (1, 1,  : ) =sqrt  (sum (squeeze (Tu)  .  '2,  1) ) ; 

%  Calcucalte  C-T 

C-T (1, 1,  : ) =ones (k_samples,  1)—1 . /squeeze (ITuI) ; 

%%  Calculate  First  order  Derivative  with  respect  to  T 
tmpl=kron_k  (u,  Tu)  ; 

two-C-T  _on_s  igma  (1,  1,  :  )  =2*  (c_T .  /  sigmasq)  ; 
tmp2  (  :  ,  1,  :  )  =kron_k  (ones  (9,1),  two-C-T_on_sigma)  ; 

A-f_dT=sum  (tmp2 .  *tmpl,  3)  ; 
clear  tmpl  tmp2 
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%%  Calculate  First  order  Derivative  with  respect  to  b 

%  Calcualte  T’Tu 

TpTu ( :  ,  1 ,  : ) =T ' * squeeze (Tu) ; 

tmp2  ( : ,  1,  :  )  =kron_k  (ones  (3,  1)  ,  two-C-T_on_s  igma)  ; 
a  -f-db=sum  (—tmp2 .  *TpTu,  3)  ; 
clear  tmpl  tmp2 

%%  First  order  Derivatice 
a  _ f-dx= [ a -f-dT  ;  A-f-db]; 

%%  Second  Order  Derivatives 
%  Hessian  Matrix  Element  H-TT 
uU-p=kron_k  (u,  U-p)  ; 
b=kron_k  ( C-T,  kron_k  ( uu.p  ,  eye  (3)  )  )  ; 

TuupTp=mult-k  (T,  mult-k  ( uu-p,  T  ' )  )  ; 
one-over -Tu-cube  (1 , 1 ,  : )  =1 .  /  (ITuI .  "3)  ; 
a=kron_k  ( one-Over-Tu-Cube  ,  kron-k  ( uu-p,  TuupTp)  )  ; 
two-on_sigma  (1,1,  : )  =2*  (1 .  /  sigmasq)  ; 

H_TT=sum  (kron-k  (two-on-sigma ,  (a+b)  ) ,  3)  ; 
clear  a  b 

%  Hessian  Matrix  Element  H_Tb 

b=kron_k  ( C-T,  kron-k  (u,  T)  +kron_k  (eye  (3)  ,  Tu)  )  ; 

a=kron_k  ( one-Over-Tu-Cube  ,  mult-k  (kron-k  (u,  Tu)  ,  mult-k  (u-p,  T  '  *T)  )  )  ; 
H_Tb=sum  (kron_k  (—two-On_sigma,  (a+b)  ) ,  3)  ; 

%  Hessian  Matrix  Element  H-bb 
TpTuupTpT=mul  t-k  (T  '  *T ,  mult-k  (uu-p,  T  '  *T)  )  ; 
b=kron_k  ( C-T ,  T  '  *T)  ; 

a=kron-k  (one-Over-Tu-Cube ,  TpTuupTpT)  ; 

H-bb=sum  (kron-k  (two-on-sigma,  (a+b)  )  ,  3)  ; 

%  Merge  Hessians 

H= [H-TT  H-Tb;  H-Tb  '  H-bb] ; 

%%  T  and  b  estimation 
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%  Calculate  new  estimate 
X_old=X; 

X=  X—H\ A-f-dx; 

%  Calculate  error  for  stop  condition 

GC_error=  (X—X_old)  '  *H*  (X—X_old) ;  %  Error  for  stop  condition  from  Alonso 

%  Extract  T=  reshape (X (1 : 9) , 3, 3) 

T=[X(1)  X(4)  X(7);  .  .  . 


X(2)  X(5)  X(8);  .  .  . 

X(3)  X(6)  X(9)  ]; 

Q, 

O 

Extract  b 

b= 

=  X  (10:12)  ; 

n= 

=n+l;  %Increase  counter 

end 

C.2.2  ‘ellipsoid  first  estimate. m’ 

function  [TO ,  bO] =ellipsoid-f irst-estimate  (h_r) 


o, 

o 

ELLIPSOID-FIRST-ESTIMATE  estimation  of  magnetometer  errors. 

o, 

o 

[TO,bO ]  =ellipsoid-f irst -estimate  (h_r) 

o 

o 

h-r  —  3xN  matrix ,  where  N  is  the  number  of  signal  samples 

g, 

o 

Noise  is  irrelevant,  and  total— field  is  assumed  to  be  constant. 

g, 

o 

g, 

o 

Implemented  based  on  the  work  presented  in: 

g, 

o 

"Geometric  approach  to  strapdown  magnetometer  calibration 

g, 

o 

in  sensor  frame " 

g, 

o 

by  J.F.  Vasconcelos  et  al .  in  Aerospace  and  Electronic  Systems, 

g, 

o 

IEEE  Transactions  47(2)  2011. 

g, 

o 

g, 

o 

This  is  a  Vasconcelos  et  al .  implementation  of: 

g, 

o 

"Extension  of  a  two-step  calibration  methodology  to  include 

g, 

o 

nonorthogonal  sensor  axes" 

g, 

o 

by  Elkaim  &  Forster  in  Aerospace  and  Electronic  Systems, 

g, 

o 

IEEE  Transactions  44(3)  2008. 

g, 

o 
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%  Note  a  correction  to  TO  ,  where  the  second  column  in  the  second  row  was 
%  chanced  to  + (1/b)  *sec  (rho)  instead  if  —  (1/b)  *sec (rho)  as  in  the  paper. 

g. 

o 

%  Justin  Dinale,  DSTO  Department  of  Defence,  Australia 
%  $Revision:  1.0.0  $  $Date:  2012/11/05$ 

0,0,0,  g,g,  0,0.0  0.0,0,0,  0.0,  o.  g,  g,  g,  g,  g,  g,  g,  g,  g,  g,  g.  g,  g,g.g.g.g.g,g,g,  g,  o_  g.  g,g.g.g.g.g,g,g,g.  0.0.  g,  g.  g.  g.  g.  g,  g,  g.  g.  Q.Q,g,g.g.g,g,g,g,g.  g.  g.  g.  g,  g. 
ooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

h_z _sq=h_r  (3,  : )  ' .  ~2; 

X=[h_r  (1,  : )  ' .  ~2,  h.r  (1,  :)  ' .  *h_r  (2,  : )  h.r  (1,  : )  ' .  *h_r  (3,  : )  ... 

h_r  (2,  :)  '  .~2  ,  h.r  (2,  : )  ' .  *h.r  (3,  : )  ... 

h.r  (1 ,  h.r  (2 ,  h.r  (3,  : )  '  ones (size  (h.r , 2) , 1) ] ; 
p=X\h.z.sq; 

A=p  (1)  ; 

B=p  (2)  ; 

C=p  (3) ; 

D=p  (4) ; 

E=p(5)  ; 

G=p  (  6)  ; 

H=p  (7)  ; 

I=P(8)  ; 

J=p  (9)  ; 

beta. 1=2 *B *H  +  B*E*I  -  2*C*D*I  -  4*D*G  +  E*C*H  -  E~2*G; 
beta.2=—2*A*E*I  +  4*A*H-  B*C*I  -  2*B*G  +  C~2*H  -  C*E*G; 
beta. 3=4  *D*I *A  -  2*D*G*C  +  E*G*B  -  I*B~2  -  2*E*H*A  +  C*B*H; 
alpha. 1  =  -B~ 2  +  D*C~2  +  4*D*A  +  A*E'~2  -  B*E*C; 
alpha.2  =  4*A*  (E~2)  *J  -  E~2*G"2  -  4*B*E*C*J  +  2*E*C*H*G  ... 

~h  2  *B  *E  *1  *G  —  4  *E  *H*A  *1  —  4*D*1*C*G  — C  2*H  2  ... 

+  4*D*A*I~2  +  2 *C*B*H*I  —  4*D*G~2  +  4*D*C~2*J  ... 

+4  *B  *H*G  —  4  *A  *H  2  —  B  2*1  2  —  4  *B  2*J  ~h  16*D*A*J / 
alpha.3  =  E~4*A-  C*B*E'~3  +  E~2*C~2*D  -  2*B~2*E~2  ... 

+  8*D*A*E~2  -  4*D*B~2  +  16*D~2*A; 

a=l/ (2 *alpha.l )  *sqrt  (—  (4 *D+E ~2)  *alpha.2 ) ; 
b=l/  (2 *alpha.l )  *sqrt  (—  (4 *A+C~2)  *alpha.2 )  ; 
c=l/ (2 *alpha.l )  *sqrt  ( (4 *D*A—B ~ 2)  *alpha.2 ) ; 

tan.rho=—0 .  5*  (2*B+E*C)  /sqrt  (alpha.l )  ; 
tan_phi=  (B*E—2*C*D)  /sqrt  (alpha.l  )  ; 
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tan_lambda=E *sqrt  (—alpha-1  / alpha-3 )  ; 
rho=atan  (tan_rho ) ; 
phl=atan  ( tan-phi) ; 
lambda=atan  ( tan_lambda) ; 

T0= [1/a,  0,  0;  ... 

—  (1/a)  *tan-rho,  +  (1/b)  *sec  (rho) ,  0;  ... 

(1/a)  *  (tan-rho  *tan-lambda  *sec  (phi)  —  tan-phl ) ,  .  .  . 

—  (1/b)  *  (sec  (rho)  *tan_lambda  *sec  (phi) ) ,  ... 

(1/c)  *  (sec  (lambda)  *sec  (phi)  )]; 

b0=  (1/ (2  *alpha-l  ))*  [beta-1  ;  beta-2;  beta-3]; 


C.2.3  ‘normcol.m’ 


function  out=normcol (A) 

%NORMCOL  Columnwise  norm. 

o, 

o 

%  normcol (A)  returns  a  vector  with  dimension  size  (A, 2)  containing 

%  the  norm  of  the  columns  of  A. the  skew  symmetric  matrix  defined  by  the 

o, 

o 

%  Example:  A=[al  a2] 

%  normcol  (A  )  =  [\\ al  \  \  ||a2||j; 

o, 

o 

%  2008—10—16  J.  F.  Vasconcelos  DSOR  SVN  repository  version. 


out=sqrt  (sum  (A.  ~2, 1) ) ; 


C.2.4  ‘mult_k.m’ 

function  M  =  mult-k(A,B) 

%MULT-K  Matrix  product ,  for  i  samples 
%  MULT-K(X,Y)  is  the  matrix  product  of  X  and  Y. 

%  The  result  is  a  large  matrix  formed  products  between  the  samples  of  X 
%  and  those  of  Y  at  each  index  i 
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o. 

o 

Q_ 

O 

■§  Help  text  adapted  from  Mathworks  KRON  function. 

%  Copyright  2012  Justin  Dinale,  DSTO  Department  of  Defence ,  Australia 
%  $Revision:  1.0.0  $  $Date:  2012/11/05$ 

9-9-9-  9-9-  9-9-9-9-9-9-9-9-9-  9-  9-9-9-  9-  9-  9-  9-  9-  9-  9-  9-9-  9-9-9-9-9-9-9-9-9-  9-9-  9-9-9-9-9-9-9-9-9-  9-9-  9-  9-  9-  9-  9-  9-  9-  9-9-  9-9-9-9-9-9-9-9-9-9-9-  9-9-  9- 
oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

sA=size (A) ;  %  Get  size  of  array  A 
sB=size (B) ;  %  Get  size  of  array  B 

sA_l=length ( sA) ;  %  Get  number  of  dimensions  for  A 
sB_l=length (sB) ;  %  Get  number  of  dimensions  for  B 

nDim=max  ( sA_l , sB_l ) ;  %  Check  what  number  of  dimension  are  used, 
if  nDim>3  %  Only  up  to  3D  are  acceptable 

error (' KRON_K  operates  with  at  most  3  dimensions'); 
elseif  nDim<3  %  Is  the  multiplication  a  single  sample? 

M=A*B; 

else  %  3D  matrix  (2D  with  multiple  samples ) 

%  Fix  up  dimension  issues  if  there  are  any 
if  sA_l<nDim 

sA= [ sA  ones (1 , nDim—sA_l ) ] ; 
k=sB  (3)  ; 

A=repmat  (A,  [1  1  k] ) ; 
elseif  sB-l<nDim 

sB=[sB  ones ( 1 , nDim—sB_l ) ] ; 
k=sA  (3)  ; 

B=repmat  (B,  [1  1  k ] ) ; 

else 

k=sA  (3)  ; 

if  (sA(3)t  sB  (3)  ) 

error (' Uneven  number  of  sampled  in  A  and  B'); 

end 

end 


%  Get  transpose  of  A. 
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A=permute  (A,  [2  1,3]); 

%Expand  out  A  for  dot  product 

A=A  ( : ,  sort  (mod  (  (1 :  sA  (1)  *sB  (2) )  —1,  sA  (1)  )  +1)  ,  : )  ; 

%  Prepare  matrix  B 

B=B  ( : ,  mod  (  (1  :  sB  (2)  *sA  (1)  )-l,  sB  (2)  )  +1,  : )  ; 

%  Perform  the  multiplication  of  A  and  B  elements . 
M=A.  *B; 

M=sum  (M,  1 ) ; 

M=reshape (M,  sB  (2) ,  [] , k) ; 

M=permute  (M,  [2  1,3]); 
end 


C.2.5  ‘kron  k.m’ 


function  K  =  kron_k_new  (A,  B) 

%KRON_K  Kronecker  tensor  product,  for  i  samples 
%  KRON_K(X,Y)  is  the  Kronecker  tensor  product  of  X  and  Y. 

%  The  result  is  a  large  matrix  formed  by  taking  all  possible 
%  products  between  the  elements  of  X  and  those  of  Y  at  each  index  k 

%  For  example,  if  X  is  2  by  3  by  k,  then  KRON_K(X,Y)  is 

o, 

o 

%  [X(1,1,:).*Y  X  (1, 2,  : )  .  *Y  X(1,3,:).*Y 

%  X  (2,  1,  : )  .  *Y  X(2,2,:).*Y  X(2,3,:).*Y] 

o, 

o 

%  No  optimisation  is  made  for  sparse  matricies  since  intended  use  is  for 
%  matricis  up  to  9x9xN 


Help  text  adapted  from  Mathworks  KRON  function. 

Copyright  2012  Justin  Dinale,  DSTO  Department  of  Defence,  Australia 
$Revision:  1.0.0  $  $Date:  2012/11/05$ 


sA=size  (A) ;  %  Get  size  of  array  A 
sB=size  (B) ;  %  Get  size  of  array  B 
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sA_l=length (sA) ;  %  Get  number  of  dimensions  for  A 
sB_l=length (sB) ;  %  Get  number  of  dimensions  for  B 

nDim=max  ( sA_l ,  sB_l ) ;  %  Check  what  number  of  dimension  are  used, 
if  nDim>3  %  Only  up  to  3D  are  acceptable 

error (' KRON_K  operates  with  at  most  3  dimensions'); 
elseif  nDim<3  %  Is  the  multiplication  a  single  sample? 

K=kron  (A, B) ; 

else  %  3D  matrix  (2D  with  multiple  samples ) 

%  Fix  up  dimension  issues  if  there  are  any 
if  sA_KnDim 

sA= [ sA  ones (1 , nDim—sA_l ) ] ; 
k=sB  (3)  ; 

A=repmat  (A,  [1  1  k] ) ; 
elseif  sB_l<nDim 

sB=[sB  ones ( 1 , nDim—sB_l ) ] ; 
k=sA  (3)  ; 

B=repmat  (B, [1  1  k] ) ; 

else 

k=sA  (3)  ; 

if  (sA(3)t  sB  (3)  ) 

error (' Uneven  number  of  sampled  In  A  and  B ' ) ; 

end 

end 

%  Prepare  matrix  A  for  multiplication .  Increase  each  element  of  A  top  be 
%  sB(l)  by  sB(2)  matrix 

A=A  (sort  (mod  (  (1 :  sA  (1)  *sB  (1)  )  —  1,  sA  (1)  )  +1) ,  .  .  . 

sort  (mod  (  ( 1 :  sA  (2)  *sB  (2)  )  —  l,  sA  (2)  )  +1) ,  : )  ; 

%  Prepare  matrix  B  using  indeximg  to  perfoprm  same  task  as  repmat 
%  B=repmat  (B,  [sA  (1) ,  sA  (2) , 1] ) ; 

B=B  (mod  (  (1  :sB  (1)  *sA  (1)  )—l,  sB  (1)  )  +l,mod  (  (1 :  sB  (2)  *sA  (2)  )—l,  sB  (2)  )  +1,  :)  ; 

%  Perform  the  multiplication  of  A  and  B  elements . 

K=A .  *B; 
end 
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C.3  AC  calibration  estimator 


C.3.1  ‘AC_estimate.m’ 


function  [D-est ,  Cov-est  ]  =AC -estimate  (B,  H,  Sigma-noise) 

%AC-ESTIMATE  estimation  of  magnetometer  errors. 

%  [D_est ,  n,  Cov_est  ]  =AC_estimate  (B,  H,  Sigma-noise  ) 

%  B  —  3xN  matrix,  where  N  is  the  number  of  signal  samples 
%  H  —  Earth  ' s  magnetic  field  in  Earth  frame  coordinates 
%  Sigma^noise  —  3xN  noise  covariance  array,  where  each  column  is  a 
%  diagonal  of  the  noise  covariance  matrix. 

g, 

o 

%  Implemented  based  on  special  case  of  the  work  presented  in: 

%  " Complete  Linear  Attitude— Independent  Magnetometer  Calibration" 

%  by  R. Alonso  and  M.  Shuster 

%  in  "The  Journal  of  the  Astronautical  Sciences"  Vol .  50  No . 4 

%  October— December  2002. 

Q, 

O 

%  Justin  Dinale,  DSTO  Department  of  Defence,  Australia 
%  $Revision:  1.0.0  $  $Date:  2013/03/19$ 

0.0.0,  0,0,  0,0,  o,  o,  q,  o,  o,  o,  g,  g,  g,  Q.Q,g,g,g,g,g,g,g,  g,  g,  g,  g,  o,Q,Q,o,g,g,g,g,g,  o,g,g,g,g,g,g,g,  g,  g,  Q,o,o,o,Q,Q,Q,o,o,o,  g,  g,  o,g,g,g,g,g,g,g,g,  g,  g,  g, 
oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo 

%%  AC  estimate 

%  Set  initial  guess  for  b  and  D. 

D0=zeros (3) ; 

%  Calculate  K 
K_k  =  [  (B .  ~2)  ;  .  .  . 

+2*B  (1 ,  :)  .  *B  (2,  :);... 

+2  *B  (1 ,  :)  .  *B  (3,  :);... 

+2  *B  (2 ,  :  )  .  *B  ( 3 ,  : )  ]  / 

%  Calculate  scalar  measurement  z 
z,k=total-field (B)  .  ~2—  total-field (H) .  ~2; 

%  Calculate  mean  and  variance  of  noise 
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mu-k^sum  (Sigma-noise,  1)  ;  %  Mean 
sigmasq-k=4  *sum  (  (  (eye  (3)  +D0)  *B)  .  * .  .  . 

(Sigma-noise .  *  (  (eye  (3)  +D0)  *B)  )  ,  1)  .  .  . 

+5 *  (sum  ( Sigma^noise  .  ~2,  1)  )—  (sum  (Sigma-noise ,  1 )  .  ~  2)  ; 

%  Calculate  fisher  information  matrix 
[ I-fisher,  I-fishinv]  =AC-fisher  ( sigma-sq-k ,  K-k)  ; 

%  Calculate  AC  estimate 
theta-0= . . . 

I-fishinv *  (—K-k  *  (  ( Z-k—mu-k )  ' .  / sigma-sq-k  ' )  )  ; 

[D_est  ]  =theta-tO-D  (theta.0 )  ; 

%  Extract  covariance  matrix 
D=D_est ; 

M-ED=  [2*D  (1) ,  0,  0,  2*D  (4) ,  2*D(5),  0;  ... 

0,  2*D  (2) ,  0,  2*D(4),  0,  2*D(6);  ... 

0,  0,  2  *D  (3)  ,  0,  2  *D  (5) ,  2*D(6);  ... 

D  (4)  ,  D  (4)  ,  0,  D  (1)  +D  (2)  ,  D(6),  D(5);  ... 

D  (5) ,  0,  D  (5)  ,  D  (6)  ,  D  (1)  +D  (3)  ,  D(4);  ... 

0,  D  (6)  ,  D  (6) ,  D  (5)  ,  D  (4) ,  D(6)+D(5)];  ... 

dD-dE=  (2*eye  ( 6)  +M-ED )  \  eye  ( 6)  ; 

Cov-est=dD-dE*  (I-fisher )  \dD-dE  ' ; 

%END  TWO STEP 

function  [B-total ] =total-f ield (B-in) 

%  Calculates  total  field  of  B-in  (3xn  matrix) . 

B_total=sum  (B_in  .  '2,  1)  .  ''0.5; 

function  [  I-fisher ,  I-fishinv]  =... 

AC-fisher  (sigma-sq,  K) 

%  Calculates  the  fisher  information  matrix  for  the  centered  estimate,  when 
%  given  variance  sigmasq  and  centered  vectors  of  K 
%  Calculate  fisher  information  matrix 
I_fisher= .  .  . 
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(  (K.  *  (ones  (size  (K,  1)  ,  1)  *  (1 .  / sigma.sq)  )  )  *K  ' )  ; 

%  Calculate  inverse 
I_fishinv=  (I_fisher)  \  eye  (6)  ; 

function  [D]=  theta.tO-D  (theta) 

%  Converts  a  value  of  theta  to  usable  physical  values 
[E]  =theta-to-E  (theta)  ; 

[U,  S]  =eig  (E)  ; 

W^eye  (3)  +  (eye  (3)  +S)  .  "(0.5); 

D=U*W*U' ; 

function  [E] =theta.tO-E  (theta) 

%  Extracts  E  elements  from  theta  as  per  Alonso  paper. 
E= [theta  (1),  theta  (4),  theta  (5);... 
theta  (4),  theta  (2),  theta  (6);... 
theta  (5),  theta  (6),  theta  (3)]; 


C.4  Parameter  recovery 


C.4.1  ‘recover  angles. m’ 


function  [ S, rho, phi,  lambda ] =recover_angles  (Dj) 

%  Convert  to  Diagonal  Matrix  for  comparison  of  estimators 
[Ql,Ll]=ql (eye  (3) +Dj) ; 

D  =  diag (sign (diag (LI) ) ) ; 

L  =  D*L1; 

Q  =  Q1*D; 

S(1)=L(1,1); 

rho=atan  (L  (2,  1)  /L  (2,  2)  )  ; 

S  (2)  =L  (2,  1)  /sin  (rho)  ; 

phi=atan  (L(3, 1)/L(3, 3)); 

lambda=atan  (L  (3,2)/  (L  (3,  3)  /cos  (phi)  )); 

S  (3)  =L  (3,  2)  /sin  (lambda)  ; 
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C.4.2  ‘ql.m’ 


function  [Q,L]=ql(A) 
tmpl=zeros  (size  (A) ) ; 
tmp2=tmpl ; 
tmpl  ( : ) =1 :numel  (A) ; 
tmp2  ( : ) =n umel (A) :  — 
J=tmpl  '  ==tmp2 ; 

[JQJ,  JLJ]  =qr  (J*A  *J)  ; 
Q=J  *JQJ  *Jf 
L=J  *JLJ 


Page  194 


Bibliography 


AAROE-M.,  MONACO-R.,  KOSHELETS-V.  R,  and  MYGIND-J.  (2009).  A  quantitative  investigation  of 
the  effect  of  a  close-fitting  superconducting  shield  on  the  coil  factor  of  a  solenoid.  Superconductor 
Science  Technology,  22(9),  p.  095017. 

ALONSO-R.,  AND  SHUSTER-M.  D.  (2002a).  Twostep:  A  fast  robust  algorithm  for  attitude-independent 
magnetometer-bias  determination.  The  Journal  of  the  astronautical  sciences,  50(4),  pp.  433-451. 
eng. 

ALONSO-R.,  AND  SHUSTER-M.  D.  (2002c).  Complete  linear  attitude-independent  magnetometer  cali¬ 
bration,  The  Journal  of  the  astronautical  sciences,  50(4),  pp.  477-490.  eng. 

ALONSO-R.,  AND  SHUSTER-M.  D.  (2003).  Centering  and  observability  in  attitude-independent 
magnetometer-bias  determination.  The  Journal  of  the  astronautical  sciences,  51(2),  pp.  133-141. 
eng. 

Auster-H.  U.,  Fornacon-K.  H.,  Georgescu-E.,  Glassmeier-K.  H.,  and  Motschmann-U. 
(2002).  Calibration  of  flux-gate  magnetometers  using  relative  motion.  Measurement  Science  and 
Technology,  13(7),  p.  1124. 

Barker-J.  R.  (1949).  New  coil  systems  for  the  production  of  uniform  magnetic  fields.  Journal  of  Scien¬ 
tific  Instruments,  26(8),  p.  273. 

Barker-R.,  Cox-M.,  FORBES-A.,  and  Harris-P.  (2007).  Software  sopport  for  metrology  best  practice 
guide  no.  4  -  Discrete  modelling  and  experimental  data  analysis.  Technical  report,  National  Physical 
Laboratory,  Teddington,  Middlesex. 

BOYVATN-M.,  AND  HAFNER-C.  V.  (2012).  Molding  the  flow  of  magnetic  field  with  metamaterials: 
magnetic  field  shielding ,  Electromagnetics  Research ,  Progress  In,  126,  pp.  303-316. 

BRACKEN-R.,  SMITH-D.,  AND  BROWN-R  (2005).  Calibrating  a  tensor  magnetic  gradiometer  using  spin 
data.  Technical  report,  U.S.  Geological  Survey,  Reston,  Virginia. 

BRONAUGH-E.  (1995).  Helmholtz  coils  for  calibration  of  probes  and  sensors:  limits  of  magnetic  field 
accuracy  and  uniformity.  Electromagnetic  Compatibility,  1995.  Symposium  Record.  1995  IEEE  In¬ 
ternational  Symposium  on,  pp.  72  -76. 

CAMPS-F.,  HARASSE-S.,  AND  MONIN-A.  (2009).  Numerical  calibration  for  3-axis  accelerometers  and 
magnetometers,  Electro/Information  Technology,  2009.  eit  '09.  IEEE  International  Conference  on, 
pp.  217 -221. 


Page  195 


Bibliography 


CAPRARI-R.  S.  (1995).  Optimal  current  loop  systems  for  producing  uniform  magnetic  fields.  Measure¬ 
ment  Science  and  Technology,  6(5),  p.  593. 

CARTER-R.  (1976).  Coil-system  design  for  production  of  uniform  magnetic  fields.  Electrical  Engineers, 
Proceedings  of  the  Institution  of,  123(11),  pp.  1279  -1283. 

Cheng-D.  (1989).  Field  and  wave  electromagnetics,  Addison-Wesley  series  in  electrical  engineering, 
Addison- Wesley. 

CLARKE-D.  (2011).  Re:  Ternan  coils.  Personal  communication. 

CLAYCOMB-J.,  AND  MlLLER-H.  (2006).  Superconducting  and  high-permeability  shields  modeled  for 
biomagnetism  and  nondestructive  testing.  Magnetics,  IEEE  Transactions  on,  42(6),  pp.  1694  -1702. 

CROTTI-G.,  AND  GlORDANO-D.  (2004).  Evaluation  of  frequency  behaviour  of  coils  for  reference  mag¬ 
netic  field  generation.  Precision  Electromagnetic  Measurements  Digest,  2004  Conference  on,  pp.  28 
-29. 

CROTTI-G.,  Chiampi-M.,  and  GlORDANO-D.  (2006).  Estimation  of  stray  parameters  of  coils  for  refer¬ 
ence  magnetic  field  generation.  Magnetics,  IEEE  Transactions  on,  42(4),  pp.  1439  -1442. 

Donley-E.  A.,  HODBY-E.,  HOLLBERG-L.,  and  KlTCHING-J.  (2007).  Demonstration  of  high- 
performance  compact  magnetic  shields  for  chip-scale  atomic  devices.  Scientific  Instruments,  Re¬ 
view  of,  78(8),  p.  7. 

DRAKE-A.  (1994).  Traceability  for  low  level  low  frequency  magnetic  field  measurements.  Low  Level 
Low  Frequency  Magnetic  Fields,  IEE  Colloquium  on,  pp.  3/1  -3/2. 

EVERETT-J.,  AND  OSEMEIKHIAN-J.  (1966).  Spherical  coils  for  uniform  magnetic  fields.  Journal  of  Scien¬ 
tific  Instruments,  43(7),  p.  470. 

Fleet  Command  -  RAN  Ranges  and  Assessing  Unit  (2009).  Commonwealth  of  Australia. 

FRIX-W.,  KARADY-G.,  AND  Venetz-B.  (1994).  Comparison  of  calibration  systems  for  magnetic  field 
measurement  equipment.  Power  Delivery,  IEEE  Transactions  on,  9(1),  pp.  100  -108. 

Griffiths-D.  (1999).  Introduction  to  electrodynamics,  Prentice  Hall. 

Hartmann-G.  K.  (1979).  Weapons  that  wait :  mine  warfare  in  the  U.S.  Navy  /  Gregory  K.  Hartmann, 
Naval  Institute  Press,  Annapolis,  Md.  :. 

Hayes-C.,  Edelstein-W.,  SCHENCK-J.,  Mueller-O.,  and  Eash-M.  (1985).  An  efficient,  highly  ho¬ 
mogeneous  radiofrequency  coil  for  whole-body  nmr  imaging  at  1.5t,  Journal  of  Magnetic  Reso¬ 
nance,  63,  p.  622  to  628. 

HUANG-L.,  AND  JlNG-W.  (2008).  Attitude-independent  geomagnetic  navigation  using  onboard  com¬ 
plete  three-axis  magnetometer  calibration.  Aerospace  Conference,  2008  IEEE,  pp.  1  -7. 


Page  196 


Bibliography 


KlM-E.,  AND  BANG-H.-C.  (2007).  Bias  estimation  of  magnetometer  using  genetic  algorithm.  Control, 
Automation  and  Systems,  2007.  ICCAS  '07.  International  Conference  on,  pp.  195  -198. 

Kim-W.-S.,  Kim-Y.  B.,  Kim-M.-S.,  Kim-K.-T.,  Chong-Y.,  Park-P.  G.,  and  Kim-Y.  G.  (2007).  Opti¬ 
mizing  the  shielding  effectiveness  of  a  shielding  cabinet  via  FEM  simulation.  Journal  of  Magnetism 
and  Magnetic  Materials. 

KlRSCHVINK-J .  L.  (1992).  Uniform  magnetic  fields  and  double-wrapped  coil  systems:  Improved  tech¬ 
niques  for  the  design  of  bioelectromagnetic  experiments.  Bioelectromagnetics,  13,  pp.  401-411. 

KLYMOVYCH-Y.,  AND  PAJUNPAA-K.  (2004).  Precise  calibration  system  for  dc  vector  magnetometers  and 
approach  to  its  accrediting.  Precision  Electromagnetic  Measurements  Digest,  2004  Conference  on, 
pp.  506  -507. 

Kuberry-R.  (1967).  Absolute  calibration  of  the  magnetic  field  generated  by  a  small  coil.  Geoscience 
Electronics,  IEEE  Transactions  on,  5(3),  pp.  94  -96. 

LASSAHN-M.,  AND  TRENKLER-G.  (1995).  Vectorial  calibration  of  3d  magnetic  field  sensor  arrays.  In¬ 
strumentation  and  Measurement,  IEEE  Transactions  on,  44(2),  pp.  360  -362. 

Lee-S.-K.,  AND  ROMALIS-M.  V.  (2008).  Calculation  of  magnetic  field  noise  from  high-permeability 
magnetic  shields  and  conducting  objects  with  simple  geometry.  Journal  of  Applied  Physics,  103(8), 
pp.  084904  -084904-10. 

Lenz-J.,  and  Edelstein-S.  (2006).  Magnetic  sensors  and  their  applications.  Sensors  Journal,  IEEE, 
6(3),  pp.  631  -649. 

LlAO-Y.  (2011).  Phase  and  frequency  estimation:  High-accuracy  and  low-complexity  techniques,  Master's  the¬ 
sis. 

MAGER-A.  (1968).  Magnetic  shielding  efficiencies  of  cylindrical  shells  with  axis  parallel  to  the  field. 
Journal  of  Applied  Physics,  39(3),  p.  1914. 

Mager-A.  (1970).  Magnetic  shields.  Magnetics,  IEEE  Transactions  on,  6(1),  pp.  67  -  75. 

MATHWORKS.  (2012).  Mathworks  Australia  -  Procrustes  Analysis  @ONLINE. 

Meeker-D.  (2010).  Finite  Element  Method  Magnetics  Version  4.2  Users  Manual@ONLINE. 

MENGCHUN-R,  HONGFENG-P.,  Shitu-L.,  Ql-Z.,  and  RuifanG-X.  (2010).  Calibration  improvement 
of  three-axis  magnetometer  in  disturbing  magnetic  circumstance  based  on  fir  digital  filter.  Signal 
Processing  (ICSP),  2010  IEEE  10th  International  Conference  on,  pp.  247  -250. 

MERAYO-J.,  BRAUER-P.,  Primdahl-F.,  J.R.-P.,  and  O.  V.-N.  (2000).  Scalar  calibration  of  vector  mag¬ 
netometers,  Measurement  Science  and  Technology,  11,  pp.  120  -132. 

MlRZAEVA-G.,  SUMMERS-T.,  AND  Betz-R.  (2012).  A  laboratory  system  to  produce  a  highly  accurate 
and  uniform  magnetic  field  for  sensor  calibration.  Industrial  Technology  (ICTT),  2012  IEEE  Interna¬ 
tional  Conference  on,  pp.  1020  -1025. 


Page  197 


Bibliography 


NlSSEN-J.,  AND  PaulsSON-L.-E.  (1996).  Influence  of  field  inhomogeneity  in  magnetic  calibration  coils. 
Instrumentation  and  Measurement ,  IEEE  Transactions  on,  45(1),  pp.  304  -306. 

Olsen-N.,  Risbo-T.,  Brauer-P.,  Merayo-J.,  Primdahl-F.,  and  Sabaka-T.  (2001).  In-Flight  Cali¬ 
bration  Methods  Used  For  The  rsted  Mission. 

PETRUCHA-V.,  AND  KASPAR-P.  (2009).  Calibration  of  a  triaxial  fluxgate  magnetometer  and  accelerom¬ 
eter  with  an  automated  non-magnetic  calibration  system.  Sensors,  2009  IEEE,  pp.  1510  -1513. 

PlERGENTILI-F.,  CANDINI-G.,  AND  ZANNONI-M.  (2011).  Design,  manufacturing,  and  test  of  a  real¬ 
time,  three-axis  magnetic  field  simulator.  Aerospace  and  Electronic  Systems,  IEEE  Transactions  on, 
47(2),  pp.  1369  -1379. 

POOR-H.  V.  (1994).  An  Introduction  to  Signal  Detection  and  Estimation,  Springer. 

SCHAAD-T.  P.  (2009).  Nano-resolution:  Oceanic,  atmospheric  and  seismic  sensors  with  parts-per-billion 
resolution.  Technical  report,  Paroscientific  Inc.,  Redmons,  Washington. 

SCHILL-R.  A.,  AND  HOFF-K.  (2001).  Characterizing  and  calibrating  a  large  helmholtz  coil  at  low  ac 
magnetic  field  levels  with  peak  magnitudes  below  the  earth's  magnetic  field.  Review  of  Scientific 
Instruments,  72(6),  pp.  2769  -2776. 

Shifrin-V.,  Alexandrov-E.,  Chikvadze-T.,  Kalabin-V.,  Yakobson-N.,  Khorev-V.,  and  Park- 
P.  (2000).  A  standard  calibration  system  for  geomagnetometers.  Precision  Electromagnetic  Mea¬ 
surements  Digest,  2000  Conference  on,  pp.  236  -237. 

SOKEN-H.,  AND  HAJIYEV-C.  (2011).  In  flight  magnetometer  calibration  via  unscented  kalman  filter. 
Recent  Advances  in  Space  Technologies  (RAST),  2011  5th  International  Conference  on,  pp.  885  - 
890. 

Stansfield-D.  (1991).  Underwater  Electro acoustin  Transducers,  Bath  University  Press  and  Institute  of 
Acoustics. 

Tashiro-K.,  Wakiwaka-H.,  Matsumura-K.,  and  Okano-K.  (2011).  Desktop  magnetic  shielding 
system  for  the  calibration  of  high-sensitivity  magnetometers.  Magnetics,  IEEE  Transactions  on, 
47(10),  pp.  4270  -4273. 

T.  Oetiker,  H.  Partl-I.  H.,  and  Schlegl-E.  (2000).  The  Not  So  Short  Introduction  to  LATEX  2e, 
Tobias  Oetiker  and  all  Contributers  to  LShort. 

Tumanski-S.  (2007).  Induction  coil  sensors  a  review.  Measurement  Science  and  Technology,  18(3), 
p.  R31. 

Vasconcelos-J.,  Elkaim-G.,  Silvestre-C.,  Oliveira-P.,  and  Cardeira-B.  (2011).  Geometric  ap¬ 
proach  to  strapdown  magnetometer  calibration  in  sensor  frame.  Aerospace  and  Electronic  Systems, 
IEEE  Transactions  on,  47(2),  pp.  1293  -1306. 


Page  198 


Appendix  C 


Bibliography 


VCELAK-J.  (2006).  Calibration  of  triaxial  fluxgate  gradiometer.  Journal  of  Applied  Physics,  99(8), 
pp.  08D913  -08D913-3. 

WANG-X.  (2008).  Automatic  and  adaptive  correction  of  diversionary  errors  in  tri-axial  magnetometer 
using  neural  networks.  Knowledge  Acquisition  and  Modeling  Workshop,  2008.  KAM  Workshop 
2008.  IEEE  International  Symposium  on,  pp.  271  -274. 

White-L.  B.  (2012).  Detection,  estimation  and  classification  2012,  Lecture  Notes,  Adelaide  University. 

Yoda-K.  (1990).  Analytical  design  method  of  self-shielded  planar  coils.  Journal  of  Applied  Physics, 
67(9),  pp.  4349  -4353. 

YOUGUANG-G.,  JlANGUO-Z.,  Zhiwei-L.,  JlNJIANG-Z.,  Lu-H.,  AND  SHUHONG-W.  (2006).  Calibration 
of  sensing  coils  of  a  three-dimensional  magnetic  property  tester.  Magnetics,  IEEE  Transactions  on, 
42(10),  pp.  3243  -3245. 

ZHANG-X.,  AND  GaO-L.  (2009).  A  novel  auto-calibration  method  of  the  vector  magnetometer.  Elec¬ 
tronic  Measurement  Instruments,  2009.  ICEMI  '09.  9th  International  Conference  on,  pp.  1-145  -1- 
150. 

Zlobin-A.,  Ambrosio-G.,  Andreev-N.,  Barzi-E.,  Chichili-D.,  Kashikhin-V.,  Limon-P., 
Terechkine-I.,  Yadav-S.,  and  Yamada-R.  (2001).  Development  of  cos-theta  nb3sn  dipole  mag¬ 
nets  for  vlhc.  Particle  Accelerator  Conference,  2001.  PAC  2001.  Proceedings  of  the  2001,  Vol.  5, 
pp.  3427 -3429  vol.5. 


Page  199 


Page  classification:  UNCLASSIFIED 


DEFENCE  SCIENCE  AND  TECHNOLOGY  ORGANISATION 
DOCUMENT  CONTROL  DATA 


1.  PRIVACY  MARKING/CAVEAT  (OF  DOCUMENT) 


2.  TITLE 

Magnetic  Test  Facility  -  Sensor  and  Coil  Calibrations 


3.  SECURITY  CLASSIFICATION  (FOR  UNCLASSIFIED  REPORTS 
THAT  ARE  LIMITED  RELEASE  USE  (L)  NEXT  TO  DOCUMENT 
CLASSIFICATION) 


Document 

Title 

Abstract 


(U) 

(U) 

(U) 


4.  AUTHOR(S) 
Justin  Peter  Dinale 


5.  CORPORATE  AUTHOR 

DSTO  Defence  Science  and  Technology  Organisation 
506  Lorimer  St 

Fishermans  Bend  Victoria  3207  Australia 


6a.  DSTO  NUMBER 
DSTO-RR-0396 


6b.  AR  NUMBER 
AR-015-690 


6c.  TYPE  OF  REPORT 
Research  Report 


7.  DOCUMENT  DATE 
August  2013 


8.  FILE  NUMBER 

9.  TASK  NUMBER 

10.  TASK  SPONSOR 

11.  NO.  OF  PAGES 

12.  NO.  OF  REFERENCES 

U-490-6-477-1 

CEI 

CEI 

221 

65 

13.  DSTO  Publications  Repository 

http:/ / dspace.dsto.defence.gov.au/dspace/ 


14.  RELEASE  AUTHORITY 
Chief,  Maritime  Division 


15.  SECONDARY  RELEASE  STATEMENT  OF  THIS  DOCUMENT 

Approved  for  public  release 


OVERSEAS  ENQUIRIES  OUTSIDE  STATED  LIMITATIONS  SHOULD  BE  REFERRED  THROUGH  DOCUMENT  EXCHANGE,  PO  BOX  1500,  EDINBURGH,  SA  5111 _ 

16.  DELIBERATE  ANNOUNCEMENT 

No  Limitations 

17.  CITATION  IN  OTHER  DOCUMENTS  Yes 

18.  DSTO  RESEARCH  LIBRARY  THESAURUS 

Ternan,  ELFcage,  TWOSTEP,  Calibration,  Magnetic,  Coil,  Estimation 

19.  ABSTRACT 

This  report  details  the  areas  in  which  signal  processing  techniques  have  enabled  DSTO  Sydney  to  enhance  its  existing  magnetic  testing 
capability,  and  present  viable  data  processing  options  for  their  implementation.  This  entails  not  only  accurate  sensor  calibration,  but 
also  calibration  of  the  excitation  coils  used  within  the  magnetic  test  system.  Reduction  of  external  noise  influences  through  the  use  of 
high  permeability  shielding  and  filtering  are  detailed,  presenting  modelling  of  coil-shield  interactions,  and  experimental  tests  for 
various  filter  characteristics.  Various  magnetic  coil  designs,  including  DSTO's  previously  undocumented  "Ternan"  coil,  are  also 
investigated,  leading  to  the  design  of  a  new  magnetic  coil,  the  "ELFcage"  coil,  being  developed.  Field  uniformity  measurements  for  the 
"Ternan"  magnetic  coil  are  also  presented. 


Page  classification:  UNCLASSIFIED 


